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I .  INTRODUCTION 


The  development  of  numerical  models  i3  guided  by  two  conflicting  goals. 
The  first  goal  is  to  include  the  most  realistic  physics  for  the  level  of 
complexity  allowed  by  computer  requirements  and  other  practical 
considerations.  The  second  goal  is  to  construct  a  model  in  which  the 
different  physical  elements  of  the  model  are  operationally  compatible  and  lead 
to  reasonable  results .  In  an  attempt  to  achieve  the  second  goal,  it  is 
sometimes  necessary  to  compromise  part  of  the  first  goal. 

The  construction  of  the  boundary-layer  package  at  Oregon  State 
University  for  the  Air  Force  Global  Spectral  Model  of  the  Air  Force  Geophysics 
Laboratory  has  been  primarily  occupied  with  development  of  improved  physical 
anH  dynamical  modelling  (first  goal)  with  only  recent  effort  directed  towards 
overall  model  compatibility  and  performance  (second  goal) .  The  previous 
contract  period  (Mahrt  et  al.,  1984)  was  devoted  almost  entirely  to 
formulatxon  of  individual  components  of  the  model  with  an  emphasis  on  physical 
consistency.  (Note:  all  citations  in  Chapter  1  appear  in  the  reference  list 
following  Chapter  2.]  This  led  to  a  rather  original  but  simple  treatment  of 
the  vegetative  canopy  (Mahrt  et  al.,  1984)  .  An  original  two-layer  model  of 
soil  hydrology  was  developed  (Mahrt  and  Pan,  1984)  which  contained  physics 
comparable  to  the  existing  models  with  many  levels  whereas  previous  two-layer 
models  were  purely  empirical.  This  development  required  careful  partitioning 
of  the  two  layers  to  control  truncation  errors  with  respect  to  natural  time 
scales  and  required  reconsideration  of  the  soil  surface  boundary  conditions. 

A  model  of  the  atmospheric  boundary  layer  was  developed  (Troen  and 
Mahrt,  1986)  which  was  sufficiently  simple  yet  allowed  growth  of  the  boundary 
layer  due  to  both  surface  heating  and  wind  shear.  The  boundary-layer  depth 
formulation  seems  to  be  rather  robust  and  has  been  recently  tested  and  adopted 
by  the  Canadian  Atmospheric  and  Environmental  Service  (AES)  and  the  Dutch 
Royal  Meteorological  Institute  (KNMI). 

Work  under  the  present  contract  has  concentrated  on  improvement  of 
specific  aspects  of  the  model  such  as  the  implied  grid-area  averaging  of  the 
surface  fluxes,  modeling  transport  induced  by  boundary-layer  clouds,  the 
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special  circumstances  of  transport  in  the  very  stable  boundary  layer,  fluxes 
over  snow  covered  surfaces  and  soil  heat  transport.  The  present  effort  has 
also  devoted  considerable  work  to  the  interaction  of  the  boundary-layer  model 
and  the  soil  hydrology  model. 

The  basic  equations  for  the  boundary-layer  package  are  summarized  in 
Chapter  II  while  the  basic  model  development  is  motivated  in  Troen  and  Mahrt 
(1986)  .  Detailed  examination  of  the  interaction  between  the  different 
submodels  of  the  boundary-layer  package  has  concentrated  on  coupling  between 
boundary-layer  development  and  heat  and  moisture  transport  in  the  soil  model. 
Boundary-layer  development  responding  to  heat  transport  from  the  surface  has 
received  considerable  attention  in  the  literature.  However  the  large  impact 
of  surface  evaporation  and  soil  moisture  transport  on  boundary-layer 
development  has  unjustifiably  received  less  attention.  The  importance  of 
coupling  between  soil  moisture  transport  and  evolution  of  the  atmospheric 
boundary  layer  in  our  model  is  illustrated  in  Chapter  III.  In  fact  the  basic 
nature  of  the  diurnal  variation  of  the  boundary  layer  changes  completely 
between  changes  of  the  three  main  stages  of  soil  drying.  The  timing  of  the 
onset  of  various  stages  of  drying  depends  on  soil  type  as  well  as  the 
potential  evaporation  imposed  by  the  atmospheric  boundary  layer. 

The  boundary-layer  package,  mainly  without  the  improvements  of  Chapters 
IV-VI,  has  been  studied  in  global  numerical  runs  within  the  Air  Force  Global 
Spectral  Model.  These  numerical  experiments  are  reported  elsewhere  (Yang  et 
al.,  1988).  Improved  formulation  of  the  surface  exchange  coefficient  is 
developed  in  Chapter  IV  of  this  report.  The  modified  exchange  coefficient 
includes  the  statistical  influence  of  subgrid  scale  variations  of  surface 
properties.  The  associated  modifications  leads  to  an  exchange  coefficient 
which  varies  more  smoothly  with  stability  at  the  transition  between  unstable 
and  stable  conditions  and  decreases  more  slowly  with  increasing  stability  in 
the  stable  case.  Implementation  of  the  modified  exchange  coefficient  into  the 
one-dimensional  boundary-layer  package  led  to  improved  performance  especially 
with  very  stable  conditions.  However  the  improvements  were  modest. 

Improved  physical  formulation  of  the  nocturnal  boundary  layer  is 
motivated,  implemented  and  tested  in  Chapter  V  of  this  report. 


A  common 


deficiency  of  boundary-layer  models  is  strong  overestimation  of  surface 
cooling  with  very  stable  conditions.  Ad  hoc  corrections  are  often  made  which 
avoid  the  real  problem;  namely,  that  standard  boundary-layer  formulations 
erroneously  "kill"  the  turbulent  mixing  too  fast  in  very  stable  conditions 
This  inadequacy  is  corrected  here  by  using  a  larger  critical  Richardson 
numoer,  using  Xondo's  formulation  for  the  eddy  diffusivity  and  applying  the 
improved  surface  exchange  coefficient  discussed  in  Chapter  IV. 

The  work  summarized  in  Chapter  VI  develops  a  formulation  for  transpcrt 
by  shallow  boundary-layer  cumulus.  This  transport  exerts  a  major  drying 
effect  on  the  boundary  layer.  Failure  to  include  the  influence  of  cumulus- 
induced  drying  in  boundary-layer  models  will  lead  to  unrealistic  moisture 
buildup.  Comparisons  are  made  between  the  formulation  developed  in  Chapter  VI 
and  the  state  of  the  art  model  of  the  European  Centre  for  Medium-Range  Weather 
Forecasts  (ECMWF)  using  actual  data.  These  comparisons  suggest  that  the 
present  development  can  lead  to  improved  modelled  transport  within  boundary 
layers  containing  shallow  cumulus.  However  the  problem  involves  interaction 
between  the  cloud-induced  mixing  and  other  aspects  of  the  boundary-layer 
model,  and,  the  problem  may  vary  considerably  between  different  types  of 
boundary-layer  clouds.  The  present  formulation  requires  more  work  and  any 
existing  evaluation  of  the  model  must  be  considered  tentative. 
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II.  THE  SOIL  AND  BOUNDARY  LAYER  MODEL 


1.  Introduction 

In  the  present  contract,  a  complicated  set  of  equations  are  utilized 
to  model  the  atmosphere,  the  soil,  and  the  vegetated  surface.  While  the 
individual  model  components  have  been  examined  previously  (Troen  and  Mahrt, 
1986;  Mahrt  and  Pan,  1984),  the  study  of  the  interaction  among  the 
components  is  an  important  goal  in  the  present  contract.  In  addition, 
parameterizations  of  shallow  cloud  convection  and  snow-cover  are 
constructed  based  on  the  existing  model  equations.  We  first  present  the 
entire  set  of  equations  and  then  present  results  of  the  current  contract. 

The  equation  set  will  be  presented  in  Section  2  and  a  brief 
description  of  the  computational  procedure  in  Section  3.  More  detailed 
explanation  of  the  equations  can  be  found  in  the  individual  chapters  of  a 
previous  report  (Mahrt  et  ai.,  1984)  .  The  equations  for  the  atmospheric 
boundary  layer  are  given  first  because  the  effect  of  the  turbulent  mixing 
is  the  goal  of  the  entire  effort  of  the  contract.  In  order  to  clo3e  the 
system  and  calculate  the  forcing  of  the  atmospheric  variables  due  to 
turbulent  mixing,  boundary  conditions  near  the  earth' s  surface  mu3t  be 
provided.  To  obtain  these  conditions,  an  atmospheric  surface  layer 
parameterization  will  be  used.  The  exchange  of  sensible  and  latent  heat 
flux  between  the  surface  layer  and  the  underlying  surface  can  only  be 
obtained  with  a  knowledge  of  the  soil  and  ocean  surface  conditions. 
Equations  for  the  atmospheric  surface  layer,  the  soil  hydrology,  and  the 
soil  thermodynamics  used  in  the  model  are  presented  following  the  boundary 
layer  equations.  The  surface  energy  balance  is  used  to  incorporate 
radiative  heating  effects  into  both  the  boundary  layer  and  the  soil  layer 
and  is  presented  last.  From  the  computational  point  of  view,  the  order  is 
actually  reversed  as  we  must  prescribe  the  external  driving  force  first. 
This  will  become  evident  in  Section  3. 
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2 .  The  Model  Equations 


This  Section  ‘ „  divided  into  four  subsections,  each  describing 
individual  aspects  of  the  planetary  boundary  layer  (PBL)  model  and  soil  model. 
Turbulent  mixing  within  the  PBL  is  described  in  §2a,  the  surface  layer  model 
of  tf  •  atmosphere  is  given  in  §2b,  the  soil  model  is  found  in  §2c,  and  the 
surface  energy  balance  calculation  is  discussed  in  §2d. 

a .  Boundary  layer  model 

The  PBL  model  i3  as  discussed  by  Troen  and  Mahrt  (1986) .  The  model 
forecasts  the  tendencies  due  to  turbulent  mixing  of  the  potential  temperature 
(9),  specific  humidity  (q) ,  and  horizontal  components  of  the  wind.  The  set  of 
prognostic  equations  is 

=  |-(k  9 

3z  \  m  oz  J 

30  _  3  f  (39  '1 

3T  -  3^%l3z  -  Yej) 

k  [*L  y  || 

Maz  Yq|) 

Here,  only  the  vertical  diffusion  terms  due  to  boundary  layer  turbulent  mixing 
are  kept  in  the  equation  to  simplify  the  presentation.  Details  of  the 
complete  equations  may  be  found  in  Troen  and  Mahrt  (1986)  . 

The  counter-gradient  correction  (y)  is  included  in  both  (2)  and  (3) 
following  Troen  and  Mahrt  (1986),  and  is  parameterized  as  follows: 
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The  counter -gradient  cor rect ions  are  evaluated  in  terms  of  the  surface  fluxes 
of  potential  temperature  and  specific  humidity,  the  boundary-layer  depth  (h) , 
the  velocity  scale  (w3j  of  the  boundary  layer  defined  as 
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and  a  constant  C,  set  to  6.5,  as  in  Troen  and  Mahrt  (1986).  In  (6),  u*  is  the 
surface  friction  velocity  and  L  is  the  Monin-Obukhov  length.  <fm  is  a  profile 
function  which  is  specified  in  (12)  below. 


The  coefficient  of  diffusivity  for  momentum  (5^)  is 
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The  eddy  diffusivity  for  heat  is  related  to  the  eddy  diffusivity  for  momentum 
in  terms  of  the  turbulent  Prandtl  number 
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where 


In  other  terms,  the  Prandtl  number  is  assumed  to  be  a  constant  and  is 
determined  as  the  value  at  the  top  of  the  surface  layer  (z3)  using  surface 
layer  similarity  theory.  As  shown  in  Eq .  9,  the  counter-gradient  term  is 
also  absorbed  in  the  Prandtl  number.  The  profile  functions  (0m,  0^)  have 
their  normal  definition  and  will  be  defined  formally  below.  The  resulting 
prediction  equation  for  potential  temperature  and  specific  humidity  will 
therefore  not  explicitly  contain  the  counter-gradient  term  and  is  actually 
identical  in  form  to  Eq.  1  (Troen  and  Mahrt,  1986)  . 


The  boundary  layer  height  is  diagnosed  as 


h  = 


Ri  r  0ov!V(h)l 
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where  Ricr  is  the  critical  Richardson  number,  9ov  is  a  reference  virtual 
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potential  temperature,  g  is  the  gravitational  acceleration,  and  8v(/i)  is  the 
virtual  potential  temperature  at  level  h.  This  approach  to  diagnosing  the  PBL 
height  also  requires  the  specification  of  a  low  level  potential  temperature 

(8ov*l -  We  define  0OV*  in  the  following  way: 


,  stable 


(w'QyOo 


9  +  C  — — unstable 

ov  w 

s 


When  the  boundary  layer  is  unstable,  the  surface  virtual  potential  temperature 
in  (11)  is  enhanced  by  an  amount  that  is  proportional  to  the  surface  sensible 
heat  flux.  The  constant  of  proportionality  is  C  and  is  chosen  as  6.5  as  in 
Eqs .  4 ,  5 ,  and  9 . 

The  nondimens ional  profile  functions  for  shear  and  thermal  gradient  are 
defined  as  follows: 
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These  forms  are  taken  from  Businger  et  al.  (1971)  and  are  functions  of  the 
height  coordinate  (z)  and  the  Monin-Obukhov  length  scale  (L) . 


The  following  variables  are  calculated  in  the  surface  layer  and  will  be 
described  in  the  next  section: 


u*,K0'j  ,  (w  q  j  ,  L 

r\  O 


These  are,  respectively,  the  friction  velocity,  the  surface  flux  of  potential 
temperature,  the  surface  flux  of  specific  humidity,  and  the  Monin-Obukhov 
length  scale. 

b.  Surface  layer  model 

The  lowest  model  layer  is  parameterized  as  a  constant  flux  surface 
layer.  The  surface  fluxes  are  parameterized  following  Louis  (1979), 
as  follows: 
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In  Eq.  14,  u*  is  the  friction  velocity,  k  is  von  K&rm&n' 3  constant,  the 
surface  wind  speed  (|V0I)  is  evaluated  at  the  lowest  model  level  and  the 
roughness  length  (zQ)  depends  on  the  type  of  surface.  The  function  (F)  and 
the  bulk  Richardson  number  for  the  surface  layer  (Rig)  are  described  below. 

In  Eqs .  15  and  16,  the  exchange  coefficient  is  a  function  of  the  surface 

wind  speed  (|VQ|),  the  height  of  the  first  model  level  (z),  the  surface 
roughness  length  ( z 0 ) ,  and  the  bulk-Richardson  number.  The  surface  air 
potential  temperature  (0Q)  and  specific  humidity  (qQ)  are  taken  at  the  first 
model  level  while  the  surface  potential  temperature  (9g)  and  specific  humidity 
(q3)  are  obtained  from  the  surface  energy  balance.  The  surface  potential 

R/C 

temperature  is  related  to  the  surface  temperature,  Ts,  by  03  -  T3(P3/pQ)  P, 
where  p3  is  taken  to  be  1000  hPa. 


The  surface  exchange  coefficient  is 
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where 
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The  constants  b  and  b'  are  specified  as  9.4  and  4.7  respectively  while  the 
coefficient  c  is  defined  as: 
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where  C*  is  7.4  (Louis,  1979) . 


The  exchange  coefficient  is  defined  in  the  same  manner  as  in  Louis 
(1979)  so  that  the  wind  speed  factor  is  absorbed  in  it.  The  bulk  Richardson 
number  for  the  surface  layer  is  defined  as: 


g  z  (9ov-esv> 

0  IV  I2 
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(20) 


where  the  subscript  v  indicates  virtual  potential  temperature.  The  bulk- 
Richardson  number  is  a  function  of  the  height  (z) ,  the  difference  between  the 
virtual  potential  temperature  of  air  at  the  first  model  level  (0OV)  and  the 
surface  virtual  potential  temperature  (9SV)  corresponding  to  the  surface 
temperature  from  the  surface  energy  balance,  and  the  air  speed  at  the  first 
model  level  (|V0|). 
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The  length  scale  for  the  surface  layer  is  the  Monin-Obuhlcov  length. 
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(21) 


The  Monin-Obukhov  length  scale  (L)  is  defined  using  surface  variables:  surface 
virtual  potential  temperature  (0SV),  friction  velocity  (u„),  and  the  virtual 
heat  flux  at  the  surface.  The  eddy  dif fusivities  are 


K  =  u*  k  z  <()  '(|-) 
m  *  TmvL 

(22) 

, -1  z. 

(23) 

\  =  u*kzVL} 

and  are  functions  of  the  friction  velocity  (u*),  height  (2)  and  the  Monin- 
Obukhov  length  scale  (L)  .  The  dimensionless  functions  (<j>m  and  were 

defined  in  Eqs .  12  and  13. 


The  only  variables  needed  to  close  the  surface  layer  model  are  T3  and  q3 
—  these  are  available  from  the  the  surface  energy  balance  calculation  (§2d) 
and  soil  model  (§2c) ,  respectively. 

c .  Soil  model 

The  3oil  model  has  been  described  previously  by  Mahrt  and  Pan  (1984)  . 

The  soil  hydrology  i3  modeled  with  a  prognostic  equation  for  0,  here  the 
volumetric  water  content: 
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(24) 
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The  coefficient  of  diffusivity  (D)  and  hydraulic  conductivity  (K)  are 
functions  cf  the  volumetric  water  content  (Mahrt  and  Pan,  1984) .  Through  the 
extremes  of  wet  and  dry  soil  conditions,  the  coefficients  0  and  K  can  vary  by 
several  orders  of  magnitude  and  therefore  can  not  be  treated  as  constant. 
Since  the  soil  model  is  a  multi-layer  model,  a  layer  integrated  form  is 

needed: 


w>  pL  +  K<e\ 

D(6)|^,  +K(0)lt 


Eq.  25  i3  valid  for  a  layer  [zf,zi+1]  =  Az^.  At  the  surface  of  the  soil,  the 
evaporation  is  called  the  direct  evaporation.  For  direct  evaporation  (E^ir^ 
at  the  air-soil  interface  (z  »  0) ,  we  have 


K(0O)  (1-Cf)  +  1(1 -of)  (26) 


where  I  is  the  infiltration  rate  and  Of  is  the  plant  shading  factor.  The 
evaporation  can  proceed  at  a  potential  rate  when  the  soil  i3  wet  (demand 
control  stage) .  When  the  soil  dries  out,  the  evaporation  (E)  can  only  proceed 
at  the  rate  the  soil  can  diffuse  water  from  the  lower  layer  (flux  control 
stage)  in  which  case 

E<Ep 
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where  Ep  is  the  potential  evaporation  rate.  The  model  also  incorporates 
transpiration  (Et)  in  the  following  manner: 


Et  =  Ep  Of  kv^T  [Azi  g(0i)] 


r 


l 


i=l 


where  kv  is  the  plant  resistance  factor  and  n  is  taken  to  be  0.5  (Pan  and 
Mahrt,  1987)  .  The  transpiration  rate  function  g(0^)  i3  defined  as: 
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(28) 


The  transpiration  limits  8ref  and  8wilt  refer'  respectively,  to  an  upper 
reference  value  and  the  plant  wilting  factor  (Mahrt  and  Pan,  1984) .  The 
canopy  evaporation  of  free  water  (Ec)  is  formulated  as 


E  = 


E  c 

P  f 


>n 


S 


(29) 


where  S,  the  saturation  water  content  for  a  canopy  surface,  is  a 


constant  chosen  to  be  2  mm.  The  canopy  water  content  (C*)  changes  as 
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(30) 


dC 

dt 


* 


a  Precip  -  E 


c 


Precipitation  increases  the  canopy  water  content  first  while  evaporation 
decreases  C*. 


Total  evaporation  is  obtained  by  adding  the  direct  soil  evaporation,  the 
transpiration  and  the  canopy  evaporation. 


E  =  E +  E  +  E  (31) 

dir  t  c  v  7 

The  total  evaporation  cannot  exceed  the  potential  evaporation  (Ep,  defined  in 
Eq.  39) .  After  obtaining  the  evaporation,  the  "surface  specific  humidity"  qs 
is  calculated  from 


This  quantity  is  that  specific  humidity  at  the  surface  which  allows  the  bulk, 
aerodynamic  relationship  to  predict  E  given  by  (31)  and  is  used  to  transmit 
information  on  the  evaporation.  Over  water,  qg  is  the  saturated  surface 
specific  humidity. 


Soil  thermodynamics  are  treated  with  a  prognostic  equation  for  soil 
temperature  (T) : 


C(0) 


dT  _  d_ 
dt  ~  dz 


Kt<6>  S’ 


(33) 


The  heat  capacity  (C)  and  the  thermal  diffusivity  (KT)  of  the  soil  are  both 
functions  of  the  soil  water  content  (0) .  While  the  heat  capacity  (C)  is 
linearly  related  to  0,  the  coefficient  of  thermal  diffusivity  (KT)  is  a  highly 
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nonlinear  func.-on  of  9  and  increases  by  several  orders  of  magnitude  from  dry 
to  wet  soil  conditions.  The  layer-integrated  from  of  (33)  is 


3T. 

AZj  C(0.) 


Kt(0) 


3T| 

3? 


(34) 

The  upper  boundary  condition  for  the  soil  thermodynamic  model  is  the  soil  heat 
flux,  G,  an  important  component  in  the  surface  energy  balance.  It  is  found  from 


z=o 


The  system  is  closed  except  for  the  potential  evaporation,  which  is  defined  in 
the  next  section. 


d.  Surface  energy  balance 

Surface  temperature  is  determined  from  the  surface  energy  balance 
method: 


(l-a)Si  +  Ll-csT4  =  G  +  H  +  LE  (36) 

S 

where  the  first  term  on  the  left-hand  side  of  Eq.  36  is  the  downward  shortwave 
radiation  (30lar  radiation) .  The  coefficient  (X  is  the  surface  albedo  and  is 
a  function  of  surface  type.  The  second  term  on  the  left-hand  side  is  the 
downward  longwave  radiation.  The  third  term  on  the  left-hand  side  is  the 
upward  longwave  radiation  and  the  coefficient  a  is  the  Stef an-Boltzmann 
constant  (equal  to  5.6696  x  10'8  W  m-2  K-4)  .  The  first  term  on  the  right-hand 
3ide  of  Eq.  36  is  the  soil  heat  flux  defined  in  Eq.  35.  The  second  term  on 

the  right-hand  side  is  the  sensible  heat  flux.  It  is  defined  as: 
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and  is  a  function  of  the  air  density  (pQ) ,  the  specific  heat  for  air 
(Cp  =  1004  J  kg"1  K*i),  the  exchange  coefficient  Eq.  17)  and  the 

difference  between  the  surface  potential  temperature  ( 0 3 )  and  the  air 
potential  temperature  at  the  first  model  level  (80) •  The  last  term  on  the 
right-hand  side  is  the  latent  heat  flux  and  is  calculated  from  Eq.  31. 

The  potential  evaporation  (Ep)  is  calculated  using  the  surface  energy 
balance  for  the  reference  state  of  an  open  water  surface: 

(l-a)Sl  +  Ll-cT'*  =  G  +  H'+LE  (38) 

where 

E  =  p  C  (q*(T' )  -  q  )  (39) 

p  ro  h  s7 

and 

H'  =  p  c  c  or  -  e  )  (40) 

ro  p  h  v  s  o' 

The  temperature  variable  ( Ts  ' )  that  appears  in  Eqs  .  38-40  is  a  fictitious 
temperature  that  the  surface  would  have  if  the  soil  is  sufficiently  wet  to 
evaporate  at  the  potential  rate.  The  variable  qs*(T'g)  in  (39)  is  the 
saturation  specific  humidity  for  this  fictitious  temperature. 

3.  Computation  Procedures 

Computationally,  the  fictitious  surface  energy  balance  for  an  open  water 
surface  is  first  U3ed  to  obtain  potential  evaporation  (Eqs.  38-40).  On  the 
left-hand  side  of  Eq.  38  are  the  downward  shortwave  and  longwave  radiative 
fluxes  and  the  upward  longwave  radiative  flux.  On  the  right-hand  side  are  the 
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soil  heat  flux,  the  sensible  heat  flux  and  the  potential  latent  heat  flux. 

The  key  quantity  to  be  determined  in  this  equation  is  the  skin  temperature 
\TS')  that  the  surface  would  achieve  if  it  was  saturated.  Eqs .  38  and  40  are 
used  to  form  a  prediction  for  Ts '  which  is  then  used  to  predict  potential 
evaporation  from  Eq.  39  (Mahrt  and  Ek,  1984) .  Both  the  soil  heat  flux  G  and 
the  exchange  coefficient  take  on  the  value  from  the  previous  time  step. 

Using  potential  evaporation  as  an  upper  limit,  the  soil  hydrology 
package  is  updated.  Eq.  26  is  used  to  obtain  direct  evaporation  from  the 
soil-atmosphere  interface.  The  terms  on  the  right-hand  side  of  Eq.  26 
represent  the  moisture  flux  at  the  surface  and  serve  to  determine  Ep  as  well 
as  the  top  boundary  condition  for  Eq .  25.  When  the  evaporative  flux  is 
greater  than  potential,  the  potential  evaporation  is  used  both  here  and  in  Eq. 
25;  otherwise,  the  calculated  flux  is  used.  Transpiration  from  plants  is 
evaluated  using  Eqs.  21-28.  When  precipitation  occurs  it  wets  the  plant 
canopy  first.  Reevaporation  occurs  at  the  rate  given  by  Eq.  29  with  the 
conservation  equation  for  canopy  water  in  Eq.  30  (Mahrt  and  Pan.  1984;  Pan  and 
Mahrt,  1987)  . 

The  soil  thermodynamic  model  (Eq.  34)  is  used  to  obtain  the  soil  heat 
flux  using  Eq.  35.  In  the  finite-difference  form  for  Eq.  35,  an  additional 
unknown  appears;  namely,  the  3kin  temperature  Tg .  Because  Ts  is  also  an 
unknown  in  the  surface  energy  balance  (Eq.  36),  the  surface  energy  balance  can 
be  solved  at  the  same  time  as  Eqs.  34-35  (Pan  and  Mahrt,  1987) .  Once  we 
obtain  Ts  the  sensible  heat  flux  is  calculated  via  Eq.  37.  When  snowcover  is 
present,  changes  are  needed  for  the  interface  and  these  changes  are  described 
in  Appendix  C  in  Chapter  V. 

Having  obtained  Ts  using  Eq.  36  and  qs  using  Eq.  32,  we  use  the  surface 
layer  paramet erization  (Eqs.  14-16)  to  obtain  the  surface  stress,  sensible 
heat  flux  and  latent  heat  flux.  Variables  used  in  Eqs.  14-16  are  further 
defined  in  Eqs.  17-20.  In  addition,  we  calculate  the  Monin-Obukhov  scale 
height  (Eq.  21)  and  the  similarity  diffusivity  profiles  and  Kh  (Eqs.  22  and 
23) .  The  non-dimensional  shear  and  thermal-gradient  are  given  in  Eqs.  12  and 
13. 


In  the  boundary  layer  model,  we  first  determine  the  height  of  the 
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Oeandary  layer  10.  .  «•  di«nsivity  eoetfioients  —  «-  suri.ce 

as.  obtained  using  Eg, .  1-9.  Einally.  me  tendencies  of  wind  velocity 
potential  tes.per.tat.  and  apecrtic  Oa»,idity  are  caloal.ted  ...  Eg. .  ■ 
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III.  INTERACTION  BETWEEN  SOIL  HYDROLOGY  AND  BOUNDARY-LAYER  DEVELOPMENT 

1 .  Introduction 

Surface  evaporation  can  substantially  reduce  surface  heating  and 
subsequent  development  of  the  daytime  boundary  layer.  As  a  result,  boundary 
layer  development  is  quite  sensitive  to  availability  of  surface  moisture  as 
previously  demonstrated  by  McCumber  and  Pielke  (1981). 

The  interaction  among  surface  evaporation,  soil  moisture  and  boundary 
layer  development  is  quite  complex  even  in  the  cloudless  case  as  noted 
schematically  in  Figure  1 .  For  example,  the  reduction  of  boundary  layer 
development  is  partially  limited  by  negative  feedbacks.  As  surface 
evaporation  moistens  the  boundary  layer,  the  potential  evaporation  normally 
decreases  which  in  turn  reduces  the  actual  evaporation.  Exceptions  include 
the  case  of  strong  downward  entrainment  of  dryer  air  where  low  humidities  are 
maintained  in  spite  of  significant  evaporation. 

On  a  longer  time  scale,  the  surface  evaporation  may  significantly  deplete 
the  soil  moisture.  This  drying  reduces  the  surface  evaporation  even  though  it 
also  acts  to  increase  the  potential  evaporation.  The  time  scale  for  this 
process  depends  on  soil  properties  as  well  as  atmospheric  conditions. 

A  suitable  set  of  observations  which  include  both  adequate  measurements 
of  soil  variables  and  atmospheric  fluxes  axe  not  available  to  study  the 
various  stages  of  drying.  In  this  paper  we  use  a  relatively  simple  model  of 
the  soil-atmosphere  system  to  identify  the  importance  of  various  interactions 
related  to  surface  evaporation.  The  results  of  this  study  or  any  modelling 
effort  will  remain  necessarily  inconclusive  until  the  required  measurements 
become  available.  Our  goal  is  to  suggest  which  interactions  are  most 
important.  Such  information  can  assist  in  the  design  of  future  observational 
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downward 


Figure  1 .  Suspected  important  interactions  between  surface  evapotranspiration 
and  boundary  layer  development  for  conditions  of  daytime  surface 
heating.  Solid  arrows  indicate  the  direction  of  those  feedbacks 
which  are  normally  positive  (leading  to  increases  of  the  recipient 
variable).  Broken  arrows  indicate  negative  feedbacks.  Two 
consecutive  negative  feedbacks  make  a  positive  one. 
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programs  as  well  as  help  identify  parts  of  the  soil-atmospheric  modelling 
which  are  most  critical. 


A  second  goal  of  this  work  is  to  provide  a  soil-atmosphere  boundary-layer 
model  wnich  is  sufficiently  simple  to  use  in  concert  with  larger  scale 
atmospheric  models.  Recent  numerical  experiments  by  Hunt  (1985)  indicate  that 
formulations  for  soil  moisture  and  surface  evaporation  presently  used  in 
general  circulation  models  have  serious  shortcomings.  The  present  formulation 
is  somewhat  more  complicated  but  physically  more  direct. 
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2. 


The  model 


The  atmospheric  boundary- layer  model  of  Troen  and  Mahrt  (1986)  is  coupled 
to  the  soil  moisture  model  of  Mahrt  and  Pan  (1984).  The  atmospheric  model 
contains  34  levels  between  the  surface  and  four  kilometers,  although  approxi¬ 
mately  the  same  results  can  be  obtained  with  as  few  as  10  levels.  The 
boundary-layer  height  in  the  model  is  determined  using  a  diagnostic  relation¬ 
ship  based  on  a  modified  bulk-Richardson  numoer  at  each  time  step.  During  the 
day,  boundary  layer  grows  in  response  to  turbulence  generated  by  surface  heat¬ 
ing.  When  the  solar  radiation  vanishes  and  if  winds  are  weak,  the  boundary 
layer  normally  collapses  to  the  first  model  layer  (-50  m). 

The  soil  model  consists  of  a  thin  upper  layer,  5  cm  thick,  which  responds 
mainly  to  diurnal  variations  and  a  thicker  lower  layer,  95  cm  thick,  which 
participates  more  in  seasonal  changes  of  soil  water  storage.  The  potential 
evaporation  is  formulated  with  a  modified  Penman  relationship  (Mahrt  and  Ek, 
1984).  The  finite  differencing  of  the  soil  model  has  been  chosen  to  minimize 
truncation  errors.  This  choice  is  based  on  comparisons  with  higher  resolution 
versions  of  the  model  up  to  100  layers.  The  truncation  errors  for  the  two- 
layer  model,  compared  to  higher  resolution  versions,  led  to  overestimation  of 
the  evaporation  of  about  10%  for  the  case  of  clay  soil  and  only  a  few  percent 
for  the  case  of  sand.  These  errors  are  small  compared  to  other  uncertainties 
such  as  treatment  of  the  soil-air  interface.  Because  an  accurate  description 
of  moisture  transport  close  to  the  soil  surface  requires  prohibitive  vertical 
resolution,  the  modelled  surface  moisture  flux  is  overestimated.  This  over¬ 
estimation  is  compensated  by  increasing  the  air-dry  values  for  the  soil  mois¬ 
ture  content  to  0.16  to  0.25  for  sand  and  clay,  respectively.  A  10-minute 
time  step  is  used  in  all  model  runs. 
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The  soil  model  of  Mahrt  and  Pan  (1984)  has  been  generalized  to  include 


soil  heat  flux  using  the  usual  thermodynamic  relationship 


where  the  volumetric  heat  capacity  C  and  the  thermal  conductivity  K  are 
formulated  as  functions  of  soil  water  content  as  in  McCumber  and  Pielke 
(1981).  A  more  detailed  discussion  of  the  soil-thermodynamic  model  is  given 
in  Appendix  A. 

The  soil  drying  period,  and  feedback  to  the  atmosphere,  usually  extends 
over  several  days  or  even  several  weeks.  Iteration  of  one-dimensional  models 
for  such  periods  leads  to  unrealistic  buildup  of  moisture  and  heat.  This 
buildup  does  not  occur  in  the  atmosphere  because  of  clear-air  radiative  cool¬ 
ing,  horizontal  advection  of  heat  and  moisture,  and  consumption  of  moisture  by 
precipitating  systems.  Such  processes  cannot  be  sensibly  formulated  within 
the  present  framework;  instead  we  specify  a  climatic  advection  or  restoring 
term  of  the  form 

(qE  "  q)/Tq 

(6e  -  0  ) /t g 

where  q  and  9  are  the  actual  values  of  specific  humidity  and  temperature,  and 
qE  and  0E  are  pseudo  equilibrium  values.  In  the  present  study  qE  is 
specified  to  be  the  initial  conditions  described  below,  while  0E  is  speci¬ 
fied  to  be  height-independent  with  a  value  of  270  K.  Heat  buildup  was 
controlled  by  specifying  a  relaxation  time  of  tg  =  10  days  while  long  term 
moisture  buildup  was  prevented  with  a  shorter  relaxation  time  of  tq  *  1 
day.  While  advection  is  pragmatically  specified  in  this  modelling  study,  it 
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is  also  thought  to  exert  a  controlling  influence  on  evaporation,  at  least  in 
some  flow  situations  (McNaughton,  1976). 


The  atmospheric  temperature  is  initialized  with  a  constant  lapse  rate 
(6  K  km-1).  The  temperature  at  the  lowest  atmospheric  model-level  is  initial¬ 
ized  at  283.6  K.  The  initial  moisture  content  of  the  atmosphere  is  specified 
to  be  3  g  kg-1  in  the  lowest  kilometer,  2  g  kg-1  between  1  and  1.2  km, 

1  g  kg-1  between  1.2  km  and  2  km,  and  0.5  g  kg-1  above  2  km.  Both  the  initial 
wind  and  the  time-independent  geostrophic  wind  are  specified  to  be  5  m  s-1 . 

I 

*  The  initial  volumetric  moisture  content  of  the  soil  is  specified  to  be  0.42,  a 

value  which  is  saturated  with  respect  to  clay  and  super-saturated  with  respect 
to  sand  leading  to  large  percolation  through  the  bottom  of  the  sand  for  the 
first  day.  The  initial  soil  temperature  is  specified  to  be  identical  to  the 
initial  value  at  the  lowest  atmospheric  level  (283.6  R). 

The  short-wave  radiative  flux  formulation  of  Holtslag  and  Van  Ulden 
(1983)  is  applied  for  45*N  starting  with  21  June.  Albedo  for  the  Earth's 
surface  is  set  at  0.25.  For  simplicity,  we  neglect  the  change  of  soil  surface 
albedo  with  soil  drying  which  can  lead  to  significant  decreases  of  potential 
evaporation  (Van  Bavel  and  Hillel,  1976).  Downward  long-wave  radiative  flux 
is  assigned  to  be  constant  corresponding  to  a  black-body  temperature  of 
270  K.  Each  numerical  experiment  is  iterated  for  21  days  in  order  to  include 
the  important  evaporation  stages. 

In  the  next  section,  four  prototype  numerical  experiments  are  iterated 
for  sand  and  clay  soil  types  and  for  geostrophic  wind  speeds  of  5  m  s-1  or 
10  m  s-1.  Diffusivity  and  conductivity  coefficients  are  specified  following 
Clapp  and  Hornberger,  1978. 
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3 .  Drying  Stages 


Radiative  fluxes,  wind  speed,  moisture  deficit,  and  atmospheric  stability 
determine  the  potential  evaporation  which  in  turn  forces  the  actual  soil 
evaporation.  When  the  soil  is  relatively  wet,  evaporation  will  be  at  the 
potential  rate  (atmospheric  demand)  as  determined  by  atmospheric  conditions. 
When  the  soil  is  sufficiently  dry,  the  rate  of  evaporation  is  controlled  by 
the  soil  moisture  gradient  in  the  upper  part  of  the  soil.  The  various  atmo¬ 
spheric  influences  on  the  potential  evaporation  interact  with  soil  moisture  in 
a  non-linear  fashion.  Some  of  the  candidate  interactions  are  noted  in  Fig.  1. 

We  first  study  the  various  stages  of  drying  occurring  during  21 -day 
iterations  by  plotting  the  solar  noon  values  of  different  variables  over  sand 
and  clay  soils  for  the  two  different  values  of  the  geostrophic  wind  speed. 
Since  vegetation,  clouds  and  precipitation  are  not  included,  extensive  drying 
and  warming  will  result. 

The  soil  drying  and  long  term  boundary- layer  changes  can  be  divided  into 
three  main  stages.  In  the  first  stage,  the  surface  evaporation  is  at  the 
potential  rate  which  decreases  slightly  with  time  (Figs.  2-3).  In  the  second 
stage,  the  actual  evaporation  decreases  rapidly  with  time  while  the  potential 
evaporation  increases  with  time.  The  second  stage  leads  to  a  near-equilibrium 
stage  (third  stage)  where  the  evaporation  and  potential  evaporation  vary 
slowly  with  time.  The  evolution  proceeds  more  rapidly  with  sandy  soil  partly 
because  sand  has  a  larger  hydraulic  diffusivity  and  conductivity  at  high 
volumetric  water  content  and  therefore  loses  more  water  to  percolation.  The 
stages  of  drying  corresponding  to  Figs.  2-3  are  similar  to  those  in  the 
modelling  study  of  Van  Bavel  and  Hillel  (1976)  except  that  they  included  the 
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dependence  of  surface  albedo  on  soil  wetness  and  neglected  adjustment  of  the 
atmosphere  to  surface  evaporation. 

a)  First  stage:  potential  evaporation 

During  the  first  stage  when  evaporation  is  at  the  potential  rate,  ootn 
tne  specific  humidity  and  the  relative  humidity  increase  with  time  (Fig.  4) 
leading  to  a  modest  decrease  of  the  potential  evaporation.  In  most  previous 
modelling  studies  of  the  drying  stages,  the  potential  evaporation  is  held 
constant.  The  slight  decrease  of  temperature  (not  shown)  and  the  correspond¬ 
ing  decrease  of  saturation  vapor  pressure  also  cause  the  potential  evaporation 
to  decrease  during  the  first  stage  of  drying.  As  a  result  of  decreased 
evaporation,  the  surface  sensible  heat  flux  to  the  atmosphere  increases  during 
this  period  (not  shown)  even  though  the  surface  temperature  decreases. 

The  boundary  layer  grows  deeper  (Fig.  5)  each  subsequent  day  partly  due 
to  the  increase  of  sensible  heat  flux.  Part  of  the  increase  is  due  to  the 
fact  that  the  boundary  layer  grows  quickly  through  the  weakly  stratified  layer 
remaining  from  the  mixed  layer  of  tne  previous  day. 

The  wind  speed  significantly  influences  the  surface  heat  budget  and 
boundary-layer  evolution  during  the  first  stage  of  drying  since  the  surface 
evaporation  is  at  the  potential  rate  which  depends  on  the  wind  speed.  On  the 
other  hand,  the  soil  type  is  of  little  importance  since  the  evaporation  is 
determined  completely  by  the  atmospheric  demand  during  this  stage  and  the 
influence  of  soil  wetness  on  albedo  is  not  included  here. 

b)  Second  stage:  rapid  decrease  of  evaporation 

However,  the  onset  of  the  second  stage  of  drying  is  determined  by  the 
soil  type.  The  clay  soil  is  able  to  meet  the  potential  evaporation  for  five 
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Figure  5.  Evolution  of  the  noontime  boundary  layer  depth. 


or  six  days  while  for  sand  the  evaporation  falls  significantly  below  the 
potential  rate  during  the  third  day.  With  the  onset  of  the  second  stage  of 
drying,  the  evaporation  depends  mainly  on  soil  type  and  is  less  dependent  upon 
wind  speed  and  other  atmospheric  properties. 

At  the  same  time,  atmospheric  conditions  change  rapidly  at  the  beginning 
of  the  second  stage  of  drying.  The  decreasing  surface  evaporation  causes  a 
sharp  increase  of  the  surface  temperature  wnicn  in  turn  increases  the  surface 
heat  flux  and  boundary-layer  growth.  Of  special  importance  is  that  the  down¬ 
ward  entrainment  of  drier  air  from  above  the  boundary  layer  can  exceed  the 
surface  evaporation  (Fig.  6a)  leading  to  divergence  of  the  upward  moisture 
flux.  This  causes  drying  of  the  boundary  layer.  As  expected,  this  net  drying 
of  the  atmosphere  occurs  first  over  sandy  soil.  Entrainment  drying  is 
encouraged  by  the  relatively  dry  air  aloft.  Entrainment  dryinq  is  thouqht  to 
be  frequently  important  in  the  evolution  of  hiqh-plains  boundary  layers  where 
air  above  the  boundary  layer  is  often  dry  (Mahrt,  1976).  In  contrast,  with 
weaker  boundary-layer  growth  and  more  humid  air  aloft,  the  entrainment  drying 
is  relatively  unimportant  (DeBruin,  1983). 

The  warming  and  drying  cause  the  relative  humidity  to  decrease  and  the 
potential  evaporation  to  increase.  However,  because  of  increasing  control  of 
evaporation  by  the  drying  soil,  the  actual  evaporation  decreases  rapidly 
during  stage  two.  The  decrease  of  surface  evaporation  during  stage  two  causes 
major  changes  in  the  development  and  structure  of  the  boundary  layer.  For 
example,  consider  the  atmosphere  profiles  at  1000  solar  time  on  day  8  (Fig. 
6b).  Over  sandy  soil,  the  surface  evaporation  is  already  quite  small  leading 
to  large  surface  heat  flux  and  vertical  profiles  of  the  heat  flux  typical  of 
the  convective  mixed  layer.  The  heat  flux  decreases  linearly  with  height 
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Figure  6.  (a)  Evolution  of  the  surface  moisture  flux  and  the  moisture  flux 

near  the  boundary- layer  top  over  sandy  soil  with  a  geostrophic  wind 
speed  of  5  ms-1.  (b)  Vertical  profiles  of  the  heat  flux  and 

moisture  flux  at  1000  hours  on  day  8  with  a  geostrophic  wind  speed 
of  5  ms“* . 


reaching  negative  values  near  the  boundary-layer  top  due  to  downward  entrain¬ 
ment  of  warmer  air. 

In  contrast,  the  surface  evaporation  over  the  clay  soil  is  still  rela¬ 
tively  large  leading  to  smaller  surface  heating  and  thinner  boundary-layer 
depth.  The  upper  two-thirds  of  the  boundary  layer  is  characterized  by  down¬ 
ward  heat  flux  associated  with  entrainment.  This  implies  that  durinq  this 
period,  the  mixing  in  the  boundary  layer  over  clay  is  driven  primarily  by  mean 
shear  whereas  mixing  in  the  boundary  layer  over  sand  is  primarily  driven  by 
convection.  This  example  shows  how  boundary-layer  development  depends  on  soil 
type  through  the  role  of  surface  evaporation. 

c)  Stage  three:  near-equilibrium 

Eventually  the  boundary  layer  approaches  an  equilibrium  state  character¬ 
ized  by  warm  and  dry  conditions.  At  noon  the  surface  evaporation  becomes 
negligible  for  sand  and  less  than  10%  of  the  potential  rate  over  the  clay 
soil.  The  boundary  layer  is  deep,  exceeding  four  kilomecers  for  sand.  At 
this  stage  of  development,  the  depth  of  real  boundary  layers  would  normally  be 
constrained  by  synoptic  or  cloud-induced  subsidence  and/or  advection  of 
smaller  boundary  layer  depth  while  some  surface  evaporation  would  be  main¬ 
tained  by  vapor  transport  in  the  soil  and  perhaps  transpiration,  all  of  which 
are  neglected  in  these  numerical  experiments. 
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4.  Transpiration 


Realistic  modelling  of  interaction  between  the  soil  and  the  atmospheric 
boundary-layer  must  include  the  influence  of  the  vegetative  canopy.  Vegeta¬ 
tion  moderates  diurnal  variations.  Furthermore,  the  difference  between  the 
three  stages  of  drying  are  not  as  distinct  since  the  vegetation  removes  water 
from  the  deeper  root  zone  which  dries  only  slowly.  Here  the  root  zone  is 
specified  to  extend  to  the  bottom  of  the  1  meter  layer.  With  deep  rooted 
plants  the  influence  of  rapid  drying  of  the  soil  surface  is  then  less 
important . 

We  formulate  the  influence  of  the  vegetative  canopy  in  the  simplest 
possible  way  which  approximates  the  most  important  aspects  of  the  canopy; 
namely,  transpiration  and  shading  of  the  soil  surface.  These  formulations  are 
detailed  in  Appendix  B. 

With  other  conditions  the  same  as  in  Section  3,  the  presence  of  the 
canopy  shading  70%  of  the  ground  extends  the  period  of  evaporation  at  the 
potential  rate  by  several  days  for  clay  (Fig.  7);  less  soil  water  is  removed 
from  near  the  soil  surface  to  meet  the  atmospheric  demand.  The  decrease  of 
evapotranspiration  during  stage  two  is  less  compared  to  the  case  of  no 
canopy.  In  other  terms,  drying  of  the  soil  surface  does  not  substantially 
reduce  the  transpiration  rate  in  that  significant  transpiration  of  deep  soil 
water  is  maintained  during  stages  2  and  3.  As  a  further  result  of  the  trans¬ 
piration,  the  boundary  layer  is  cooler,  more  moist  and  not  as  deep  compared  to 
the  case  with  no  vegetation.  A  fourth  stage  where  transpiration  decreases  due 
to  depletion  of  deep  soil  moisture  was  not  studied  here. 
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Figure  7.  Evolution  of  noontime  evaporation  for  the  case  of  surface 
transpiration  as  compared  to  the  standard  case  with  no 
transpiration . 
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5.  Influence  of  Solar  Radiation,  Climatic  Advection,  and  Soil  Properties 

It  is  instructive  to  study  the  sensitivity  of  the  above  conclusions  to 
variations  of  the  external  forcing.  This  is  most  simply  carried  out  by 
neglecting  the  canopy.  We  first  examine  the  winter  case  (solstice)  where  the 
incoming  solar  radiation  is  much  reduced.  Under  such  conditions,  the  drying 
stages  evolve  more  slowly  (Fig.  8a)  due  to  much  lower  rates  of  potential 
evaporation.  The  second  stage  of  drying  with  sand  does  not  begin  until  after 
one  week,  while  the  evaporation  from  the  clay  soil  remains  near  the  potential 
rate  during  the  entire  twenty  one  day  period  of  numerical  integration.  The 
potential  evaporation  reaches  only  about  100  W  m~2  during  mid-day  so  that 
transport  of  moisture  within  the  clay  soil  is  able  to  meet  the  demand.  During 
the  night,  vertical  transport  within  the  clay  is  able  to  restore  the  soil 
moisture  near  the  surface  to  the  extent  that  the  evaporation  is  near  potential 
during  the  subsequent  daytime  period. 

In  actual  atmospheric  conditions,  advection  of  heat  and  moisture  can 
significantly  alter  the  boundary-layer  evolution  even  on  short  time  scales. 
Here  we  study  the  influence  of  advection  as  formulated  in  Section  2  for  the 
summertime  case.  With  less  dry-air  advection,  the  boundary  layer  moistens 
which  reduces  the  potential  evaporation  and  significantly  delays  the  transi¬ 
tion  to  the  second  stage  of  drying  for  clay,  as  is  evident  in  Figure  8b,  for 
the  case  where  the  relaxation  time  for  moisture  is  increased  to  10  days.  With 
reduced  cold-air  advection  (not  shown),  the  boundary  layer  heats  up  and  grows 
faster  which  in  turn  increases  the  downward  flux  of  dry  air. 

When  the  soil  is  thinner,  it  stores  less  moisture.  As  a  result,  the 
second  stage  of  drying  begins  slightly  earlier.  As  an  example,  decreasing  the 
soil  depth  from  1  m  to  1/2  m  advances  onset  of  the  second  stage  of  drying  by 
only  a  day  or  less,  depending  on  soil  type  (Fig.  8c).  However,  the  influence 


Latent  heat  flux  (W  m 


of  thinner  soil  becomes  more  significant  at  later  times  when  the  moisture 
content  of  the  thin  soil  decreases  to  near  air-dry  values.  The  surface 
evaporation  rates  at  stages  2  and  3  are  only  a  fraction  of  the  corresponding 
values  for  a  1  m  thick  soil. 

Often  natural  surfaces  are  covered  by  organic  "debris"  or  "litter" 
consisting  of  dead  grass  and  leaves,  conifer  needles  and  other  organic 
matter.  Such  materials  cover  a  major  portion  of  natural  land  surfaces.  When 
dry,  these  materials  are  characterized  by  extremely  low  hydraulic  conductivi¬ 
ties;  such  surfaces  then  act  as  a  moisture  barrier  and  the  soil  becomes 
decoupled  from  the  atmosphere  on  short  time  scales  in  cases  with  little  trans¬ 
piration. 

The  thermal  and  hydraulic  properties  of  such  organic  debris  are  not  known 
with  quantitative  accuracy.  Better  known  are  the  properties  of  peat  soils 
which  are  between  those  properties  of  organic  surface  material  and  those  of 
other  soils.  In  this  study,  we  use  the  thermal  and  hydraulic  properties  for 
peat  adopted  by  McCumber  and  Pielke  (1981).  The  saturation  water  content  for 
peat  is  nearly  twice  that  of  the  other  soil  types.  The  hydraulic  conductivity 
coefficient  of  saturated  peat  is  similar  to  that  of  other  soil  types  at  or 
near  saturation.  However,  as  soil  moisture  decreases,  the  hydraulic 
conductivity  for  peat  decreases  rapidly  to  values  several  orders  of  magnitude 
smaller  at  water  contents  comparable  to  the  saturation  values  of  sand  and 
clay.  At  this  stage  the  peat  becomes  an  effective  moisture  barrier  which 
decouples  direct  exchange  between  the  soil  and  atmosphere. 

Several  numerical  iterations  were  performed  where  the  upper  5  cm  was 
specified  to  be  peat  (not  shown).  The  contribution  of  the  surface  evaporation 
to  the  surface  energy  balance  quickly  becomes  negligible.  In  any  event,  the 
usual  neglect  of  organic  litter  in  large  scale  modelling  studies  probably 
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leads  to  significant  overestimation  of  evaporation  from  the  soil  over 
vegetated  natural  surfaces.  Results  are  only  useful  qualitatively  since  the 
interface  between  the  organic  material  and  the  more  conventional  soil  cannot 


be  modelled  with  certainty.  Furthermore,  organic  litter  can  reduce  run-off  by 
absorbing  more  rain  water.  This  can  actually  lead  to  increased  evaporation  at 
a  later  stage  during  near  drought  conditions. 
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6. 


Diurnal  Variation 


As  an  example,  the  diurnal  variation  of  the  surface  energy  budget  is 
shown  in  Pig.  9  for  the  standard  case  with  sandy  soil  during  days  3  and  4, 
which  corresponds  to  the  beginning  of  the  second  stage  of  soil  drying.  Note 
than  uii  uay  4,  tile  evaporation  in  significantly  reduced  leading  to  areater 
sensible  heat  flux  to  the  atmosphere. 

The  surface  evaporation  increases  during  the  morning  as  dictated  by 
increasing  net  radiation  and  resulting  increase  of  potential  evaporation 
(Fig.  9).  This  rapid  increase  of  evaporation  suppresses  sensible  heat  flux  to 
the  atmosphere  and  leads  to  temporary  retardation  of  heat  flux  to  the  soil. 

By  late  morning,  the  soil  surface  layer  has  dried  to  the  extent  that  the 
evaporation  becomes  subpotential  and  decreases  in  an  absolute  sense.  With 
less  evaporation,  the  sensible  heat  flux  to  the  atmosphere  increases  rapidly 
and  more  heat  is  transported  into  the  soil.  Note  that  the  heat  flux  to  the 
soil  peaks  before  the  sensible  heat  flux  to  the  atmosphere,  since  the  soil 
warms  more  rapidly  than  the  atmosphere.  The  delay  of  the  diurnal  increase  of 
sensible  heat  flux  to  the  atmosphere  is  often  observed  (e.g.,  Oke,  1978). 

7 .  Conclusions 

In  the  above  modelling  study,  the  soil  drying  advances  in  three  stages  as 
has  been  previously  observed.  In  the  first  stage,  the  rate  of  surface 
evaporation  proceeds  at  the  potential  rate  and  therefore  depends  mainly  on 
atmospheric  conditions  such  as  wind  speed,  relative  humidity  and  incoming 
solar  radiation.  Surface  heating  is  limited  by  the  surface  evaporation  and 
the  boundary  layer  may  develop  primarily  due  to  shear.  In  such  cases,  weak 
downward  heat  flux  can  extend  downward  throughout  much  of  the  boundary  layer. 
In  the  second  stage,  the  evaporation  decreases  rapidly  to  well-below  potential 
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values  and  becomes  controlled  more  by  the  moisture  gradients  in  the  soil.  In 
the  final  stage,  the  drying  reaches  a  small  near-equilibrium  value.  The 
surface  heat  flux  becomes  much  larger  than  the  latent  heat  flux  and  upward 
heat  flux  extends  upward  to  the  entrainment  region  of  the  boundary  layer  top. 

The  duration  of  each  stage  depends  critically  on  the  soil  type  as  well  as 
on  atmospheric  conditions.  The  occurrence  of  dry  organic  debris,  such  as 
leaves  and  dead  grass,  appears  to  partially  decouple  the  atmosphere  and  soil 
resulting  in  significant  slowing  of  the  advance  of  the  second  and  third  stages 
of  drying.  The  development  of  significant  transpiration  reduces  the 
importance  of  direct  surface  evaporation  from  the  soil  and  thus  reduces  the 
distinction  between  the  three  stages.  Entrainment  drying  of  the  boundary 
layer  can  become  important  with  dry  air  aloft  and  strong  surface  heating.  The 
latter  is  encouraged  by  dry  soil  conditions. 

Ciimatic  cooling  and  drying  was  specified  to  simulate  advection,  clear 
air  radiative  cooling  and  removal  of  moisture  by  convective  clouds.  Such 
mechanisms  are  not  necessary  when  simulating  only  one  diurnal  cycle  or  when 
the  one-dimensional  model  is  combined  with  a  larger  scale  model  We  are 
presently  using  the  boundary  layer-soil  model  with  the  global  spectral  model 
of  Brenner  et  al . ,  1984,  as  reported  further  by  Yang  et  al^. ,  1988. 
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APPENDIX  A:  The  Two-Layer  Soil  Thermodynamic  Model 


1 


•\ 


The  two-layer  structure  used  for  the  soil  moisture  model  (Mahrt  and  Pan, 
1984)  should  adequately  resolve  the  diurnal  variation  of  the  soil  thermo¬ 
dynamics;  the  thin  top  layer  with  a  thickness  of  5  cm  can  provide  an  estimate 
of  the  sharp  diurnal  thermal  gradient  and  the  thicker  second  layer  (95  cm) 
allows  us  to  incorporate  heat  storage  and  seasonal  variations  and  to  specify  a 
constant  lower  boundary  soil  temperature  which,  in  reality,  varies  on  the 
annual  time-scale. 

The  heat  conduction  equation,  neglecting  horizontal  interactions,  is 
given  as 
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3T 
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77 


IK  |Z) 

3z 
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where  C  is  the  volumetric  heat  capacity  and  K  is  the  thermal  conductivity. 
The  heat  capacity  for  water  is  4.2  x  106  W  m-3  K"1  and  for  soil  is  chosen  as 
1.26  x  106  w  m~3  K-1  for  simplicity  even  though  it  varies  slightly  for 
different  soil  *-ypes.  The  heat  capacity  of  the  composite  soil  is  simply 
defined  as 


C  =  (1 


0)  C 


soil 


+  ec 


water 


where  0  is  the  volumetric  water  content.  In  this  definition,  we  have 
neglected  the  contribution  due  to  air  following  DeVries  (1975).  The  thermal 
conductivity,  K,  is  strongly  dependent  on  the  soil  moisture  content.  Similar 
to  McCumber  and  Pielke  (1981),  we  also  adopted  the  functional  form  for  K 
following  A1  Nakshabandi  and  Kohnke  (1965): 

420  exp  (-P.  +  2.7)  P  5.1 

K ( 0 )  =  {  (A-2 ) 

0.1722  P  >  5.1 


46 


where  Pf  =  log^g  [i|/s(0s/9  )b] .  The  factors  tj)s,  9S,  and  b  are  functions  of  the 
soil  textural  class  (Clapp  and  Hornberger,  1978). 

In  the  finite  difference  formulation,  the  model  equation  (A-1)  will  be 
integrated  first  over  the  two  layers  to  explicitly  express  the  flux  K  3T/3z, 
through  each  layer.  The  model  grid  staggering  is  presented  in  Fig.  A1  and  the 
level  zk  represent  the  level  along  which  the  temperature  Tk  is  the  same  as 
the  layer-average  T  (in  this  study,  the  mid-point  of  the  layer  is  selected). 
The  layer  integrated  equation  becomes 
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where  the  gradient  3T/3z  is  evaluated  as 
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At  the  top  of  the  model,  the  surface  temperature  Ts  will  be  used  to  form  a 
one-sided  estimate  of  the  gradient 


3T 

Tz 

The  gradient  at  the  bottom  of  the  model  is  estimated  using  a  specified 
constant  temperature,  Tbot  (Fig.  A1 ) . 

In  order  to  interface  the  soil  thermodynamics  into  the  model,  the 
prediction  of  Tk  using  (A-3)  is  performed  using  the  fully  implicit 
Cranck-Nicholson  scheme  given  by 
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where  the  superscripts  designate  the  time  levels. 
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Fig.  A-l.  The  geometry  of  the  soil  thermodynamics. 


Figure  A-1.  The  geometry  of  the  soil  thermodynamics. 
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For  moist  soil,  a  small  difference  in  the  thermal  gradient  in  (A-5)  can 
lead  to  large  soil  heat  flux  because  the  thermal  conductivity  increases 
rapidly  with  soil  moisture  content.  For  this  reason,  the  surface  energy 
balance  equation  must  be  updated  simultaneously  with  the  soil  thermodynamic 
equations  so  that  the  resulting  surface  temperature  and  soil  temperature 
satisfy  the  surface  energy  balance  constraint  at  the  next  time  step. 


I 
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APPENDIX  B:  Transpiration 


This  appendix  describes  the  transpiration  formulation  for  the  results  in 
Section  4.  We  want  to  preserve  the  distinction  since  the  direct  soil 
evaporation  is  most  appropriately  related  to  the  soil  moisture  of  an  upper 
thin  layer  while  water  for  transpiration  originates  more  from  the  deeper  root 
zone.  The  total  evaporation  can  be  written  as 
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E_.  +  Em  + 
dir  T 


(B-1  ) 


where  E^j.  is  the  direct  evaporation  from  the  soil,  ET  is  the 
transpiration,  and  E^  is  the  evaporation  of  precipitation  intercepted  by  the 
canopy.  Each  of  the  evaporation  terms  on  the  right-hand  side  are  proportional 
to  the  potential  evaporation  Ep  (we  do  not  differentiate  between  ground 
temperature  and  "leaf"  tenperature  in  this  study). 

Vegetation  reduces  the  direct  evaporation  from  the  soil  by  shading  the 
ground  and  reducing  the  wind  speed  near  the  ground.  The  reduction  of  wind 
speed  can  be  posed  in  terms  of  increased  surface  roughness  parameter  and 
increased  displacement  height. 

The  reduction  of  solar  radiation  reaching  the  ground  surface  through  the 
vegetation  can  be  expressed  as  a  linear  dependence  on  the  shading  factor  by 
neglecting  complexities  due  to  varying  sun  angle. 

To  minimize  the  number  of  parameters,  we  relate  both  the  influences  of 
shading  and  wind  speed  reduction  to  the  shading  factor  of  according  to  the 
format 


E_.  =  E 
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ESoii  is  the  evaporation  from  the  soil  in  the  absence  of  vegetation  as 

\ 

discussed  in  Section  2. 
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>  Transpiration  is  related  to  the  density  of  vegetation  and  the  soil 
moisture  content.  For  the  two-layer  model,  these  influences  are  most  simply 
included  with  the  following  formulation  for  transpiration 

ET  =  Epkv°f  g(9l)  +  9(92)l1  -  <C*/S)n] 
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where  Zj  is  the  depth  of  the  upper  layer  (here  5  cm)  and  z2  is  the  depth  of 
the  entire  two  layers  (1m).  We  have  assumed  that  the  root  uptake  rate  is 
independent  of  depth  within  a  given  layer.  After  consulting  numerous  studies, 
the  wilting  point,  where  root  uptake  ceases,  is  assigned  to  be  0.12. 

The  parameter  9ref  is  the  soil  moisture  content  where  the  soil  moisture 
deficit  begins  to  reduce  root  uptake  and  transpiration.  9ref  is  chosen  to 
be  0.25  which  is  significantly  below  the  saturation  values  for  most  soil 
types. 

C*  is  the  canopy  water  content  and  S  is  the  canopy  water  capacity  and  are 
included  to  represent  reduction  of  evaporation  of  heat  surfaces  covered  by  a 
water  film.  The  coefficient  kv  is  the  plant  resistance  factor  chosen  to  be 
1.0  and  Of  is  specified  to  be  0.7  The  product  of  kv<jf  is  similar  to 
the  commonly  used  plant  coefficient.  The  parameter  n  is  chosen  to  be  1/2  to 
be  consistent  with  the  interception  model  discussed  below. 

Some  dewfall  occurs  in  the  iterations  reported  in  Section  4. 

Interception  is  modelled  as 
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( B-4  ) 
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3E~  ~ 
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Ec  =  af(CVS)nEp 

The  storage  capacity  of  the  canopy,  S,  is  chosen  to  be  2  mm.  P  is  the 
precipitation  or  dewfall  rate.  This  interception  model  is  similar  to  that  of 
Rutter  et  al.  (1971)  except  that: 

1 )  the  throughf all  parameter  is  replaced  with  the  closely  related 
expression  1  -  Of  in  order  to  reduce  the  number  of  parameters; 

2)  the  evaporation  factor  C*/S  is  multiplied  by  of  to  account  for  the 
asymptotic  limit  that  canopy  evaporation  vanishes  as  the  canopy 
vanishes ;  and 

3)  n  is  chosen  to  be  less  than  unity  to  correspond  to  a  finite  time  for 
the  canopy  to  dry  following  rainfall  as  modelled  in  Deardorff  (1978). 

Based  on  the  work  of  Leyton  et  al.  (1967),  a  value  of  n  =  1/2  is  inferred 
which  is  somewhat  less  than  the  n  =  2/3  value  chosen  by  Deardorff  (1978). 

Once  the  canopy  is  saturated  (C*  =  S) ,  all  additional  rainfall  is  assumed 
to  fall  through  to  the  ground.  This  is  analogous  to  assuming  that  drip 
processes  occur  instantaneously  so  that  the  canopy  is  never  temporarily 
"supersaturated. " 
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IV.  GRID-AVERAGED  SURFACE  FLUXES 


1 .  Introduction 

Numerical  models  of  geophysical  flows  require  parameterization  of  the 
transport  by  all  motions  which  are  not  resolved  by  the  grid.  The  parameter¬ 
ization  of  such  "subgrid-scale  flux"  at  the  surface  is  normally  based  on 
boundary  layer  similarity  theory  and  definition  of  a  surface  exchange 
coefficient. 

Formulations  of  sub-grid  scale  flux  suffer  several  major  problems: 

a)  They  do  not  explicitly  include  transport  by  motions  which  are  larger  than 
turbulent  scale  but  still  small  enough  to  be  subgrid  scale.  These 
motions  include  nonlinear  gravity  waves,  cloud-induced  motions  and  flow 
responding  i-u  auLgj.id  terrain  ana  differential  heating.  Transport  by 
such  motions  is  poorly  understood  because  they  are  usually  observed  with 
significant  sampling  problems.  Since  the  smallest  resolved  motion  and 
the  largest  subgrid-scale  motions  are  of  comparable  scale  and  may  be 
strongly  interactive,  the  transport  by  the  largest  subgrid-scale  motions 
cannot  be  simply  related  to  the  resolved  gradient. 

b)  The  surface  is  inhomogeneous  on  subgrid  scales.  Because  the  transport  by 
the  turbulence  is  related  to  gradients  and  stability  in  a  nonlinear  way, 
the  area-averaged  flux  is  not  related  to  the  area-averaged  gradient  in  a 
simple  manner.  For  example  the  vertical  gradient  of  the  area-averaged 
potential  temperature  often  corresponds  to  stable  stratification  even 
though  the  area-averaged  heat  flux  is  upward;  that  is,  strong  turbulence 
in  small  regions  of  unstable  stratification  can  dominate  the 
area-averaged  heat  flux  which  then  becomes  "counter”  to  the  area-averaaed 
vertical  gradient  of  potential  temperature.  Small  subregions  of  unstable 
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stratification  could  result  from  local  terrain  elevation,  dryer  soil,  or 
cloud-free  pockets.  As  a  result  of  similar  averaging  effects,  the 
downward  grid-averaged  heat  flux  in  the  stable  case  may  be  substantially 
larger  than  predicted  by  the  stability  based  on  grid-averaged  variables. 
c)  With  stronger  subgrid-scale  inhomogeneity,  the  turbulence  may  not  achieve 
equilibrium  with  the  local  surface  in  which  case  practical  representa¬ 
tions  of  turbulence  are  not  applicable. 

The  above  averaging  problems  are  rarely  formally  recognized  in  modeling 
studies.  Wyngaard  (1982)  and  others  have  examined  the  mathematics  of  grid- 
volume  averaging  for  cases  where  the  grid  volume  is  both  larger  and  smaller 
than  the  characteristic  scale  of  the  turbulence.  Sud  and  Smith  (1984) 
simulate  idealized  grid-volume  averaging  by  assuming  that  the  surface  bulk 
Richardson  number  varies  within  a  grid  box  according  to  a  Gaussian  frequency 
distribution.  The  surface  exchange  coefficient  for  the  bulk  aerodynamic 
relationship,  which  depends  on  the  Richardson  number,  is  then  averaged  over 
the  hypothetical  grid  volume.  The  resulting  averaged  exchange  coefficient 
exhibits  a  much  smoother  dependence  on  stability  which  eliminates  troublesome 
numerical  oscillations.  Further  application  of  the  smoothed  exchange 
coefficient  is  found  in  Sud  and  Smith  (1985). 

Analogous  problems  have  been  studied  with  respect  to  longer  term  time 
averaging.  Saltzman  and  Ashe  (1976a, b)  have  considered  contributions  to  the 
monthly  averaged  heat  flux  due  to  diurnal  and  synoptic  variations  and  how  such 
contributions  relate  to  a  local  flux-gradient  formulation.  Mahrt  et  al. 

(1986)  have  studied  certain  averaging  problems  associated  with  use  of  the 
surface  flux  relationships  with  long  time  steps  or  omission  of  the  diurnal 
variation.  For  the  data  sets  examined,  the  actual  long-term  surface  heat  flux 
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was  upward  and  counter  to  the  long-term  vertical  gradient  of  potential 
temperature. 

In  the  present  study,  we  examine  the  influence  of  horizontal  averaging  on 
the  relationship  between  the  flux  and  gradients.  In  particular  we  horizontal¬ 
ly  average  the  local  surface  flux  relationship  to  show  how  the  nonlinear 
dependence  of  the  exchange  coefficient  on  local  gradients  generally  leads  to 
larger  flux  in  the  stable  case  than  would  be  predicted  by  the  usual  neglect  of 
spatial  averaging. 

The  formal  grid-area  averaging  of  the  flux-gradient  relationship  is 
developed  in  Section  2.  In  Section  3,  the  surface  exchange  coefficient  is 
averaged  for  idealized  distributions  of  the  Richardson  number.  Low-level 
aircraft  observations  in  the  stable  boundary  layer  are  considered  in  Section 
4.  In  Section  5,  various  averaging  terms  are  evaluated  by  using  a  mesoscale 
model  and  viewing  it  as  a  grid  box  of  a  larger  scale  general  circulation 
model.  A  modified  relationship  for  the  surface  exchange  coefficient  is 
constructed  in  Section  6. 


2. 


Formulation 


In  numerical  models,  the  flow  is  automatically  divided  into  the  resolved 
part  and  the  unresolved  or  subgrid  part.  The  flux  divergence  due  to  the 
subgrid  flow  influences  the  resolved  flow  and  must  be  parameterized.  In  this 
study  we  are  concerned  with  fluxes  from  the  earth's  surface  through  the  lowest 
atmospheric  level  in  the  model. 

We  partition  the  flow  at  the  lowest  model  level  as 

$(x,y,t)  =  [$1  +  $(x,y,t)  (la) 

where  [<H  is  the  grid  area-average  of  the  local  time  averaged  part  of  the  flow 

[$]  =  7T  /  /t+T  $(x,y,t)  dt  dA  ,  (1b) 

A  dA  t 

a 

and  is  the  deviation  from  the  grid-averaged  part.  Here  the  independent 
variables  x  and  y  refer  to  position  within  the  grid  area,  and  t  refers  to  a 
"fast"  time  for  averaging  the  local  flow.  This  averaging  eliminates  turbulent 
fluctuations  corresponding  to  time  scales  smaller  than  x.  The  grid-averaged 
flow  is  constant  in  terms  of  these  independent  variables,  but  of  course  varies 
on  larger  space  and  time  scales. 

In  order  to  apply  existing  formulations  for  the  surface  fluxes,  the 
subgrid-scale  flow  must  be  partitioned  into  a  local  time-averaged  part  £*(x,y) 
and  fluctuations  from  this  time  average  $'(x,y,t)  so  that 

*(x,y,t)  =  $*(x,y)  +  $'(x,y,t)  (1c) 

where  .  , 

1  t-Hr 

$*(x,y)  =  —  /  $(x,y,t)  dt  -  [$]  .  (Id) 

t 

The  part  $ '  is  usually  referred  to  as  turbulent  fluctuations  although  in 
practice  $ '  includes  all  motions  whose  time  scales  are  smaller  them  the 
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averaging  time  t.  Existing  formulations  for  the  surface  fluxes  are  based  on 
fluxes  estimated  from  observations  which  use  averaging  times  typically  between 
ten  and  thirty  minutes.  The  flow  component  $>*  includes  all  motions  on  spatial 
scales  larger  than  the  "turbulent"  scale  (time  scale  t)  and  smaller  than  the 
scale  of  the  resolved  flow.  In  global  models  with  grid  resolutions  of  100  Jem 
or  more,  includes  mesoscale  motions.  The  particular  partitioning  of  the 
flow  selected  above  is  not  the  only  possibility  but  provides  a  useful 
framework  for  identifying  the  important  grid-averaging  problems. 

The  expression  for  the  grid  averaged  surface  flux  using  the  above  decom¬ 
positions  becomes 

[wifi]  =  [  w  ]  [<t>  ]  +  [w*<|>*]  +  [w'4>’]  +  [  [w]if>*] 

+  [w*[<H]  +  [  [  w]<t> '  ]  +  [  w  ’  [  4>  ]  ]  +  [w*<j>  '  ]  +  [  w '  4>  *  ] 

The  last  six  terms  on  the  right-hand  side  are  cross  terms  which  vanish  for 
simple  unweighted  time  and  space  averaging;  that  is,  [w]  can  be  pulled  outside 
the  space  and  time  integrals  defining  the  operator  (  ] ,  w*  can  be  pulled 
outside  the  time  integral,  the  areal  integral  cf  w*  vanishes  and  the  time 
integral  of  w’  vanishes.  With  weighted  averaging  such  as  filtering,  the  cross 
terms  would  not  normally  vanish  (e.g.,  Charnock,  1957).  With  simple 
unweighted  averaging,  the  expression  for  the  grid-averaged  surface  flux 
simplifies  to 

[w$  ]  =  [w][<j>]  +  [w*$*]  +  [w '  4>  '  ]  (2) 

The  term  [ w ] [ <^  ]  is  the  resolved  flux  which  is  usually  converted  to  advection 
format  by  applying  incompressible  mass  continuity  to  the  divergence  of  the 
flux.  The  term  [w*$*]  represents  the  vertical  flux  due  to  the  time-averaged 
flow  which  varies  spatially  within  the  grid  area.  Near  the  surface,  w*  is 
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normally  small  except  in  terrain- related  flows.  Higher  in  the  boundary  layer, 
the  transport  term  [w*<j>*]  can  be  large  partly  because  the  spatial  (wavenumber) 
energy  gap  between  the  turbulent  scales  included  in  $ '  and  the  resolved  flow, 
[w],  often  disappears.  The  subgrid-scale  flux  [w*$*]  is  so  situation 
dependent,  there  is  no  practical  way  to  parameterize  it  in  terms  of  the 
resolved  flow.  As  a  result,  sophisticated  formulations  for  the  remaining  flux 
are  not  justified  for  use  in  large-scale  models. 

To  apply  existing  relationships  for  the  turbulent  surface  fluxes  in 
familiar  form,  we  symbolize  the  time-averaging  operator  with  an  overbar  and 
using  (Id)  define  the  total,  local,  time-averaged  flow  as 

1  t+T 

<Mx,y)  =  —  /  <Mx,y,t)  dt  =  [  ]  +  4>Mx,y) 

T  t 

It  will  also  be  useful  to  express  [w'$ 1 ]  as 

—  J  /t+T  w 1  $  '  dt  dA  =  -[  /  w 1 1 j>  1  dA  =  [  w '  <f> '  ]  . 

TA  dA  t  A  dA 

The  use  of  the  overbar  is  redundant  but  poses  the  grid-averaging  problem  in 
terms  of  the  usual  local  operators. 

With  these  expressions,  we  can  now  relate  the  surface  turbulent  flux  to 
the  local  time-averaged  flow  using  the  usual  bulk  aerodynamic  relationships 
for  surface  fluxes.  The  basic  problem  for  numerical  models  is  that  such 
existing  formulations  relate  the  local  turbulent  flux  to  the  local  vertical 
gradient.  In  numerical  models  it  has  been  necessary  to  use  the  identical 
relationship  for  relating  the  grid-averaged  flux  to  the  vertical  gradient  of 
the  grid-averaged  (resolved)  flow  even  though  there  is  no  justification  for 
such  use.  Except  for  the  exploratory  study  of  Sud  and  Smith  (1984),  there  are 
no  existing  formulations  for  relating  the  area-averaged  flux  to  the  flow 
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gradients  resolved  by  the  model.  The  following  analyses  will  show  that  the 
difference  between  the  local  flux-gradient  relationship  and  the  grid-averaged 
one  is  more  than  a  mathematical  subtlety. 

Although  it  is  not  possible  to  construct  a  relationship  for  [ w*<f>  *  J ,  it  is 
possible  to  examine  plausible  behavior  of  the  area  average  of  the  flux  due  to 

the  "turbulence"  [  w  ’  tf>  '  ] .  For  example,  consider  the  local  surface  flux 
formulated  with  the  bulk  aerodynamic  relationship,  which  in  present  notation 
is  of  the  form 

iwT)  ,  =r  v(i  -  $  > 
sfc  $  sfc 

or,  in  terms  of  the  model  resolved  flow 

(wT)  -  =  ( [c.  ]  +  C*)((V]  +  V*){([$  ,]+♦*,)-  (t<fr]  +  <*.*)}  (3) 

sfc  <p  <t>  sfc  sfc  1 

where  "sfc"  refers  to  the  surface  value  while  other  variables  are  defined  with 
respect  to  the  lowest  atmospheric  level,  is  the  surface  exchange 
coefficient,  V  is  the  surface  wind  speed,  and  $  represents  the  transported 
quantity  such  as  heat,  moisture  or  momentum.  The  surface  exchange  coefficient 
depends  on  the  choice  of  the  atmospheric  level,  such  as  10  m  or  50  m,  and  is 
especially  sensitive  to  the  stability  of  the  flow.  The  validity  of  this 
relationship  (3)  in  evolving  boundary  layers  is  discussed  in  the  next  section. 

Spatially  averaging  the  surface  flux  relationship  over  the  grid  volume, 
we  obtain 

[(wT^sfci  =  (vmaw  -  iti)  +  [c;v*]([*sfc]  -  un 

+  [C  )IV*(**fc  -  ♦*>]  +  [V]  [C*(4>*fc  -  ♦*)!  +  [C*  V*{|*fc  -  **}]  (4) 

where,  aaain,  only  area-averaged  variables  are  resolved  by  numerical  models. 
Substitution  of  this  expression  into  (2)  then  defines  the  total  vertical  flux 


of  <p  within  a  grid  box.  However,  only  the  first  term  in  (2)  can  be  computed 
in  modeling  studies  and  then  only  with  the  additional  approximation  that  the 
spatially  averaged  exchange  coefficient  is  related  to  that  stability  evaluated 
from  spatially-averaged  variables.  Formally, 

1  =  ^  -  Sj,(S'  z0>  dA  (5) 

whereas  numerical  models  can  evaluate  only 

vif  [Zol>  (6) 

where  S  is  the  stability  parameter  and  Zq  is  the  surface  roughness  parameter; 
dependencies  on  the  boundary-layer  depth  and  thermal  wind  are  not  considered 
here.  The  tilde  signifies  that  the  function  is  computed  from  variables  that 
are  already  averaged  over  the  entire  grid  box  as  opposed  to  averaging  the 
function  itself.  Sud  and  Smith  (1984)  have  examined  the  behavior  of  the  area 
average  of  the  surface  exchange  coefficient  of  Deardorff  (1972)  for  a  special 
case  where  the  wind  speed  and  roughness  was  constant  over  the  grid  area  and 
the  subgrid  variation  of  vertical  temperature  gradient  obeyed  a  Gaussian 
distribution.  In  their  case,  [C^ ]  =  .  Except  for  their  study,  the 

dependence  of  [C^ ]  on  stability  has  received  little  formal  attention. 

In  summary,  the  formulation  of  subgrid  fluxes  in  a  numerical  model  suffer 
three  types  of  errors:  a)  omission  of  the  subgrid  flux  [w*$*],  b)  omission  of 
the  various  interaction  terms  in  the  expression  for  the  averaged  surface  flux 
(4),  and  c)  approximation  of  [C^ ]  in  terms  of  area-averaged  variables 
instead  of  area-averaging  the  exchange  coefficient.  Analogous  averaging 
problems  occur  with  formulation  of  fluxes  at  model  levels  above  the  surface. 

In  the  naxt  section,  we  examine  the  spatial  averaging  of  the  surface  exchange 
coefficient. 
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3. 


Area-Averaged  Exchange  Coefficient 


The  area-averaged  exchange  coefficient  may  be  quite  different  from  the 
exchange  coefficient  computed  from  the  stability  parameter  based  on  area- 
averaged  variables.  This  is  because  the  turbulence  and  exchange  coefficient 
depend  on  the  stability  in  a  nonlinear  way.  As  one  example,  consider  the  case 
where  the  area-averaged  stratification  is  stable  but  varies  within  the  grid 
area.  The  exchange  coefficient  predicted  by  the  area-averaged  stability  may 
be  quite  small  in  stable  conditions.  However  due  to  the  strong  nonlinearity 
of  the  stability  dependence,  the  area-average  of  the  local  exchange 
coefficient  may  be  significantly  larger  due  to  subgrid  areas  where  the 
stratification  is  near-neutral  or  unstable.  In  these  subgrid  areas,  the  local 
exchange  coefficient  may  be  one  or  more  orders  of  magnitude  larger  tnan 
implied  by  tne  spatially  averaged  stability.  Because  of  the  nonlinear 
dependence  of  the  excnange  coefficient  on  stability,  small  subgrid  regions 
could  have  a  strong  influence  on  the  grid-averaged  exchange  coefficient  and 
f.lux  but  little  influence  on  the  grid-averaged  stability. 

(a)  Relationship  for  exchange  coefficient 

To  illustrate  such  averaging  problems,  we  will  adopt  the  formulation  of 
Louis  (1979)  for  the  exchange  coefficient  for  heat.  This  formulation  closely 
approximates  similarity  theory  but  is  considerably  simpler  and,  consequently, 
has  found  considerable  use  in  large  scale  models.  The  exact  form  of  the 
relationship  for  the  exchange  coefficient  is  not  too  important  provided  that 
it  includes  the  rapid  decrease  of  the  exchange  coefficient  with  increasing 
stabi’ity  which  occurs  at  near-neutral  values  of  stability.  The  behavior  of 
the  exchange  coefficient  at  strong  stability  is  quite  uncertain  partly  due  to 
flux  sampling  problems  such  as  those  discussed  by  Wyngaard  (1973).  The 
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ability  to  assign  roughness  values  and  the  occurrence  of  nonequilibrium 
conditions  over  realistic  surfaces  also  reduces  the  importance  of  the  details 
of  the  exchange  relationship. 

The  Louis  formulation  relates  the  surface  exchange  coefficient  to  the 
surface  bulk  Richardson  number1 

Ri  =  [0 ( z )  -61  z/(V(z))2  (7) 

sfc 

where  9sfC  is  the  potential  temperature  corresponding  to  the  surface 
temperature,  z  is  the  neight  of  the  first  model  level  above  ground  and  0Q  is  a 
basic  state  temperature  scale.  The  dependence  of  the  Louis  excnange 
coefficient  is  plotted  in  Figure  1  for  three  different  values  of  z/zg  .  The 
rapid  variation  of  the  coefficient  at  near-neutral  stability  will  dominate  the 
averaging  effects. 

(b)  Spatial  variation  and  adjustment 

The  surface  exchange  coefficient  varies  primarily  due  to  interrelated 
spatial  changes  of  surface  temperature,  wind  speed  and  surface  roughness. 

Many  land  surfaces  experience  continuous  changes  of  surface  characteristics 
whereas  previous  studies  have  concentrated  primarily  on  discontinuities  of 
surface  properties.  Previous  studies  have  also  largely  neglected  variations 
of  surface  evapotranspiration  which  can  occur  over  surfaces  which  otherwise 
appear  to  be  homogeneous.  This  inhomogeniety  is  forced  by  variations  of 
vegetation  and  variations  of  soil  moisture.  Soil  moisture  variations  are 
forced  by  spatial  changes  of  soil  type  and  small  scale  precipitation 

1  Although  the  Richardson  number  approaches  infinity  with  free  convection,  the 
dependence  of  the  Louis  exchange  coefficient  on  the  Richardson  number  is 
such  that  the  predicted  heat  flux  remains  well  behaved  in  the  free 
convection  limit. 
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patterns.  Such  inhomogeneity  can  force  significant  spatial  variations  of 
surface  heat  flux,  atmospheric  stability  and  boundary  layer  development 
(Ookouchi  et  _al . ,  1984;  Pan  and  Mahrt,  1987).  As  a  result,  all  surfaces  are 
potentially  inhomogeneous  and  the  problem  of  adjustment  to  surface  inhomogene¬ 
ity  must  be  considered  before  applying  existing  flux-gradient  relationships. 

Existing  formulations  for  surface  fluxes  invoke  a  form  of  the 
flux-gradient  relationship;  nonequilibrium  formulations  do  not  exist.  Even 
though  the  fluxes  and  other  boundary  layer  characteristics  may  change  rapidly 
downstream  from  a  change  of  surface  properties,  existing  models  assume  that  a 
local  equilibrium  is  maintained  betweeen  the  local  surface  flux  and  the  local 
vertical  gradient.  For  example  in  studies  of  neutral  flow  over  a 
discontinuity  of  surface  roughness,  the  flux-gradient  relationship  is  assumed 
in  the  form  of  the  logarithmic  relationship  as  the  lower  boundary  condition 
(Peterson,  1969;  Panchev  et_  al.  ,  1971;  Rao  et_  al . ,  1974;  and  others). 

Observational  studies  of  flow  over  changes  of  surface  heating  and 
roughness  have  normally  concentrated  on  the  horizontal  variation  of  the  depth 
of  the  internal  boundary  layer.  Little  attention  has  been  devoted  to  the 
degree  of  internal  equilibrium  between  surface  fluxes  and  local  vertical 
gradients.  Some  studies  of  internal  boundary  layers  generated  by  onshore  flow 
have  shown  that  turbulence  statistics  at  relatively  short  distances  from  the 
shore  show  considerable  agreement  with  similarity  relationships  based  on 
observations  over  homogeneous  terrain.  Such  similarity  describes  much  of  the 
data  collected  by  Smedman  and  Hbgstrbm  (1983)  1500  m  from  the  shore  in  weakly 
heated  onshore  flow.  In  a  study  of  sea  breeze  flow  2  km  inland,  Mizuno  (1982) 
also  found  turbulence  statistics  to  be  described  by  local  similarity  theory, 
although  the  depth  of  the  internal  boundary  layer  exerted  a  greater  influence 
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compared  to  observations  over  homogeneous  terrain.  The  local  equilibrium 
implied  by  the  success  of  similarity  theory  was  argued  by  Mizuno  in  terms  of 
small  values  of  the  turbulent  adjustment  time  scale  compared  to  the 
Langrangian  time  scale  of  the  mean  flow  between  the  coast  and  observation 
point . 

For  example  the  adjustment  length  scale  for  turbulence  equilibrium  can  be 
posed  as 

L  =  U  l/u 

where  l  is  the  length  scale  of  the  main  eddies,  probably  some  fraction  of  the 

depth  of  the  internal  boundary  layer,  u  is  the  turbulence  velocity  scale  wnich 

can  be  estimated  as  the  square  root  of  the  turoulence  kinetic  energy  and  U  is 
a  scale  vaiue  for  the  speed  of  the  mean  flow.  For  plausible  values  of 

4=100  m,  u=1  m/s,  and  U=5  m/s,  the  adjustment  length  scale  is  500  m.  For  this 

particular  example,  the  transition  region  where  equilibrium  conditions  are  not 
approximately  valid,  would  be  narrow  compared  to  the  grid  width  of  most  larger 
scale  models.  However  it  is  not  known  how  to  apply  this  information  to 
surfaces  with  continuous  changes  of  surface  conditions.  Surface  features 
which  are  smaller  scale  than  the  main  tranporting  motions  in  the  boundary 
layer  will  probably  not  exert  an  important  influence  on  the  overall  boundary 
layer  flux.  The  influence  of  such  surface  variations  will  be  integrated  by 
the  main  boundary-layer  eddies. 

While  it  is  not  possible  to  consider  realistic  surfaces  where  distinct 
internal  boundary  layers  are  often  prevented  by  complex  continuous  changes  of 
surface  conditions,  it  is  necessary  that  the  idealized  subgrid  variations  are 
consistent  with  the  application  of  the  local  flux-gradient  approximation. 

More  specifically  it  must  be  assumed  that  changes  in  surface  conditions  are 
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either  sufficiently  gradual  that  the  flux-gradient  relationship  remains  a 

>  useful  approximation  or  that  nonequi librium  regions  occupy  a  small  fraction  of 
the  grid  area.  Then  the  bulk  aerodynamic  relationship  is  a  useful 
approximation  if  one  can  account  for  subgrid  variations  of  the  exchange 
coefficient  due  to  horizontal  variations  of  stability,  roughness,  and  boundary 
layer  depth.  Here  we  address  the  variations  of  the  exchange  coefficient  due 

► 

to  variations  of  stability  which  appears  to  be  the  most  important  influence. 

i 

f  (c)  Distribution  and  averaging 

Variations  within  a  given  grid  area  depend  on  geograpnic  location,  time 

of  day,  season,  and  synoptic  situation.  It  is  not  possible  nor  practical  to 
i 

t 

explicitly  consider  such  variations  in  large  scale  models.  Here  we  consider 
the  probability  distribution  of  a  generic  grid  box.  A  Gaussian  probability 
density  of  the  Richardson  number  is  suitable  for  demonstrating  the  potential 
importance  of  averaging  errors. 

Employing  the  Gaussian  distribution  as  in  Sud  and  Smith  (1985),  the 
averaged  exchange  coefficient  becomes 

[C]  =  1  f  C  (Ri)dA  =  /  C  (Ri)  f ( Ri ) dRi  (8) 

where  f(Ri)  is  the  assumed  Gaussian  probability  density  of  the  Richardson 
number  over  area  A  where  Ri  (7)  is  defined  in  terms  of  local  time-averaged 
variables.  Relationship  (8)  assumes  nothing  about  the  spatial  coherence  or 
pattern  of  the  Richardson  number  but  assumes  that  the  overall  distribution  is 
Gaussian.  The  integral  on  the  right-hand-side  of  (8)  was  evaluated  using 
Simpson's  rule.  It  was  found  that  using  a  step  size  of  .  1o  over  a  range  of 
±8o  provided  accurate  results  in  that  additional  widening  of  the  range  or 
shortening  the  step  size  no  longer  significantly  altered  the  results. 
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As  an  example,  Figure  1  shows  the  area-averaged  surface  exchange 
coefficient  for  a  Gaussian  distribution  of  the  Richardson  number  with  unity 
standard  deviation  and  zero  mean.  The  variation  of  the  Richardson  number 
together  with  the  nonlinear  dependence  of  the  exchange  coefficient  cause  the 
spatially  averaged  exchange  coefficient  to  be  significantly  larger  for  stable 
conditions  compared  to  the  values  for  the  original  local  relationship.  In 
other  words,  when  the  area-averaged  Richardson  number  is  stable,  the  spatially 
averaged  exchange  coefficient  may  be  dominated  by  the  small  part  of  the  area 
corresponding  to  the  unstable  tail  of  the  frequency  distribution.  In  the 
sub-areas  of  near  neutral  or  unstable  stratification,  the  exchange  coefficient 
is  much  larger  than  in  stable  areas  and  therefore  has  an  important  influence 
on  the  area-averaged  value  even  if  most  of  the  area  is  stable. 

The  influence  of  averaging  is  minimal  for  near-neutral  average  stability 
where  the  dependence  ' i  the  exchange  coefficient  on  the  Richardson  number  is 
characterized  by  an  inflection  point.  The  averaging  influence  on  the  exchange 
coefficient  is  percentage-wise  small  for  large  instability  where  the 
dependence  of  the  exchange  coefficient  on  the  Richardson  number  becomes  more 
linear.  Considering  the  uncertainty  of  the  original  formulation  for  the 
exchange  coefficient,  the  influence  of  averaging  for  the  above  conditions  is 
probably  important  only  for  the  stable  case. 

The  enhancement  of  the  exchange  coefficient  for  stable  conditions  is  even 
greater  with  larger  standard  deviation  of  the  Richardson  number  since  larger 
standard  deviation  implies  that  a  greater  portion  of  the  area  will  be  near 
neutral  or  unstable.  This  is  indicated  in  Figure  2,  where  the  area-averaged 
exchange  coefficient  is  shown  for  different  values  of  the  standard  deviation. 

When  the  standard  deviation  of  the  Richardson  number  exceeds  about  1/2, 
the  averaging  effect  begins  to  completely  change  the  behavior  of  the  exchange 
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Figure  1.  Dependence  of  the  averaged  exchange  coefficient  on  the  averaged 

Richardson  number  for  three  different  values  of  z/z0  for  standard 
deviation  of  the  Richardson  number  equal  to  one  (thick  lines)  and 
the  original  unaveraged  exchange  coefficient  (thin  lines). 


Figure 


.  Dependence  of  the  averaged  excnan< 
deviation  of  the  Richardson  numbe] 
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coefficient.  This  change  is  partly  due  to  extreme  values  at  deviations 
greater  than  one  or  two  standard  deviations  toward  the  unstable  regime.  When 
the  standard  deviation  of  the  Richardson  number  becomes  very  large,  say 
greater  than  five,  the  dependence  of  the  averaged  exchange  coefficient  on  the 
averaged  Richardson  number  becomes  almost  linear  and  the  rapid  change  of  the 
exchange  coefficient  at  near-neutral  stability  disappears. 

The  averaging  effect  could  be  even  greater  with  significant  skewness  of 
the  distribution  of  the  Richardson  number.  However,  typical  distributions  of 
the  Richardson  number  over  a  given  land  area  are  not  known.  In  the  next 
section,  limited  information  on  frequency  distributions  of  the  Richardson 
number  is  computed  from  data  collected  with  low-flying  aircraft. 
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4. 


Observations 


We  have  computed  the  Richardson  number  from  data  collected  with 
low-flying  aircraft  in  the  stable  nocturnal  boundary  layer  over  gently 
undulating  terrain  in  south  central  Oklahoma,  USA,  during  quiet  periods  of  the 
Severe  Environmental  Storms  and  Mesoscale  Experiment  (SESAME) .  The  aircraft 
instrumentation  is  described  in  Mahrt  (1985)  and  Wyngaard  et  al.  (1978). 

The  Richardson  number  is  computed  between  the  aircraft  level  and  the 
surface.  The  surface  temperature  is  estimated  from  the  radiation  temperature 
inferred  from  the  downward  pointing  radiometer  mounted  on  the  aircraft.  The 
estimated  error  of  the  surface  temperature  is  thought  to  be  small  compared  to 
the  temperature  difference  except  tor  near-neutral  conditions.  Errors  due  to 
the  assumption  of  unity  surrace  emissivity  are  partially  cancelled  by  the 
surface  reflection  of  downward  longwave  radiation.  More  importantly,  the 
difference  between  the  surface  land  temperature  and  the  air  temperature  at  z  = 
zq  must  be  neglected,  as  in  most  modeling  situations.  Such  differences  may  be 
several  degrees  with  strong  cooling  or  heating. 

The  surface-based  layer  Richardson  number  in  analogy  with  (7)  is 


Ri 


(g/eG  ><e-9rad>z 
V2 


(9) 


where  the  potential  temperature  9  and  wind  speed  V  are  averaged  over  segments 
of  the  aircraft  record  and  0racj  is  the  potential  temperature  corresponding 
to  the  surface  radiation  temperature.  The  level  z  of  the  horizontal  flights 
ranges  between  20  m  and  100  m  for  the  various  legs. 
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To  estimate  the  "local"  average  $  in  lieu  of  a  time  average,  the  record 
is  divided  into  75  m  segments  (20  observational  points)  and  the  Richardson 
number  was  computed  from  variables  averaged  for  each  segment.  This  estimate 
is  useful  if  the  grid  width  of  the  intended  model  is  large  compared  to  75  m. 
With  increased  averaging  length,  the  standard  deviation  of  this  Richardson 
numoer  for  a  given  flight  leg  decreases  but  the  relationship  between  the 
Richardson  number,  its  standard  deviation  and  the  turbulence  intensity  does 
not  change  appreciably  for  the  data  analyzed  here. 

Statistics  for  the  spatial  distribution  of  the  Richardson  numoer  were 
computed  for  each  of  tne  37  aircraft  legs  which  were  typically  15-30  tan  long. 
The  standard  deviation  of  the  Richardson  number  within  a  given  flight  leg 
increases  significantly  with  increasing  stability  of  the  leg  (Figure  3a).  As 
a  result  of  this  variability,  some  subregions  of  weak  stability  and 
significant  turbulence  are  possible  even  when  the  averaged  stability  is 
large.  As  a  further  consequence,  the  turbulence  and  turbulence  flux  are 
unlikely  to  vanish  even  with  large  averaged  stability.  Note  that  in  Figure 
3a,  the  standard  deviation  is  plotted  as  a  function  of  the  mean  of  the 
Richardson  numbers  computed  along  the  flight  leg.  This  mean  Richardson  number 
was  generally  closely  related  to  the  mean  Richardson  number  computed  from 
variables  averaged  along  the  leg,  although  one  can  imagine  realistic 
situations  where  this  relationship  would  not  hold. 

The  dependence  of  the  standard  deviation  of  the  Richardson  numoer  on  its 
mean  value  does  not  seem  to  be  sensitive  to  the  aircraft  flight  level  although 
both  the  standard  deviation  and  mean  value  tend  to  decrease  with  the  height  of 
the  aircraft  level.  This  decrease  is  due  to  the  decrease  of  stratification 
with  heif.it.  The  correlation  between  the  mean  value  and  standard  deviation  of 
the  Richardson  number  is  partly  due  to  the  fact  that  the  Richardson  number  is, 
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with  a  few  exceptions,  bounded  by  zero  since  the  surface  is  almost  everywhere 
cooler  than  the  overlying  air  for  this  data.  As  a  result,  greater  standard 
deviation  of  the  Richardson  number  leads  to  greater  mean  value  and  vice  versa. 

The  Richardson  number  computed  from  the  present  data  is  often  character¬ 
ized  by  skewness  towards  large  positive  values  although  this  behavior  is  too 
erratic  to  incorporate  into  the  analysis  of  Section  3.  On  one  of  the  days, 
the  wind  speed  nearly  vanishes  leading  to  extremely  large  positive  Richardson 
numbers  which  in  turn  cause  the  standard  deviation  and  mean  value  of  the 
Richardson  numbers  to  be  "off  scale"  for  Figure  3.  This  occurred  in  three  of 
the  37  aircraft  legs.  This  behavior  is  one  of  the  natural,  but  unfortunate, 
ciiaracteristics  cf  the  Richardson  number  which  can  lead  to  misleading 
statistics. 

Since  the  standard  deviation  of  the  Richardson  number  increases  with 
stability,  the  most  appropriate  averaging  for  the  results  reported  in  Figure  2 
would  use  larger  standard  deviations  at  larger  Richardson  numbers.  This  has 
the  effect  of  slowing  the  decrease  of  the  exchange  coefficient  with  increas¬ 
ing  stability  as  composited  over  many  distributions  with  different  means.  The 
data  cannot  be  used  to  directly  study  the  relationship  between  the  actual  flux 
and  the  layer  Richardson  number.  Fluxes  computed  at  the  flight  level  would  be 
contaminated  by  large  sampling  problems.  Since  the  fluxes  are  relatively  weak 
and  intermittent  under  very  stable  conditions,  a  much  longer  record  over 
relatively  homogeneous  terrain  would  be  needed. 

However,  it  was  possible  to  compute  the  variance  of  the  vertical  velocity 
as  an  indicator  of  turbulence  strength  for  each  record  segment.  The  variance 
is  less  sensitive  to  sampling  problems.  The  variance  computed  for  each  75  m 
segment  is  undoubtedly  due  almost  exclusively  to  small  scale  turbulence  with 
little  influence  of  gravity  waves.  The  variances  were  then  further  averaged 
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(a)  Standard  deviation  of  the  Richardson  number  for  each  aircraft 
leg  as  a  function  of  the  leg-averaged  Richardson  number  for  the 
evening  of  4  May  1979  (open  circles),  early  morning  of  5  May 
(squares),  evening  of  5  May  (triangles),  early  morning  of  6  May 
(solid  circles),  early  morning  of  7  May  (crosses),  ar.d  early 
morning  of  9  May  (solid  squares).  (b)  Composited  vertical 
velocity  variance  as  a  function  of  the  Richardson  number. 


for  different  stabilitycl  isses  basea  on  the  value  of  the  Richardson  number; 
each  class  corresponds  to  a  Richardson  number  interval  of  0.25. 

The  dependence  of  the  vertical  velocity  variance  on  the  stability 
suggests  three  distinct  regimes  (Fig.  3b).  For  weak  stability  (Ri  <  1),  the 
strength  of  the  turbulence  does  not  appear  to  be  sensitive  to  the  value  of  the 
Richardson  numoer  although  the  class,  0  <  Ri  <  0.25,  contains  fewer  cases  and 
may  be  subject  to  inadequate  sampling.  With  moderate  stability  (1  <  Ri  <  3 ) , 
the  turbulence  decreases  linearly  with  increasing  stability.  In  the  very 
stable  case  (Ri  >  3),  the  turbulence  is  quite  weak  and  not  very  sensitive  to 
the  strength  of  the  stability.  Apparently,  when  the  turbulence  is 
sufficiently  suppressed  by  the  stratification,  some  residual  weak  turbulence 
remains  regardless  of  the  srreng t.\  of  the  stability.  Presumably,  this  weak 
turbulence  occurs  or.  scales  smaller  than  that  used  to  comDute  the  Richardson 
number  although  the  nature  u£  motion  was  not  studied. 

The  critical  values  of  the  Richardson  number  corresponding  to  the 
transitions,  and  the  existence  of  the  sharp  transitions  themselves  may  depend 
upon  the  way  m  which  the  Richardson  numer  is  defined.  It  must  be  noted  that 
the  Richardson  number  used  here  is  a  layer  Richardson  number  in  contrast  to 
the  local  gradient  Richardson  number  wnere  the  turbulence  is  thought  to  be 
suppressed  for  values  greater  than  about  0.25. 

For  the  present  data,  the  transition  values  did  not  depend  on  the  choice 
of  averaging  length  used  to  define  the  local  average.  Using  only  the  lowest 
aircraft  levels  (20-35  m),  the  transitions  are  sharper  and  shifted  toward 
slightly  smaller  values  of  the  Richardson  number.  With  higher  flight  levels, 
the  transitions  are  less  defined  and  shifted  to  significantly  larger  values  of 
the  Richardson  number.  That  is,  with  stable  stratification,  the  turbulence  at 
higher  levels  becomes  more  determined  by  conditions  at  that  level  and  less 
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related  to  surface  processes.  Furthermore,  Richardson  numbers  computed  over 
thicker  layers  may  be  large  but  still  allow  for  turbulence  in  thinner 
sublayers . 

Although  the  relationship  between  the  turbulence  variance  and  the 
stability  depends  on  the  way  in  which  the  Richardson  numoer  is  computed,  the 
above  turbulence-stability  regimes  are  similar  to  those  found  by  Hondo  et  al. 
(1978,  Fig.  6).  The  main  difference  is  that  in  the  study  ol  Hondo  et  al. 
(1978)  the  transition  to  the  very  stable  regime  of  weak  turbulence  occurs  at 
smaller  Richardson  numbers.  This  is  probably  related  to  the  fact  that  in 
their  study,  the  Richardson  number  was  computed  over  thinner  layers. 

In  conclusion,  the  surface-based  Richardson  number  is  a  useful  indicator 
of  turbulence  strength  in  spite  of  the  fact  that  the  Richardson  number  varies 
dramatically  at  low  wind  speeds.  More  sophisticated  formulations  of  the 
averaging  problem  (8)  might  take  advantage  of  the  well  defined  relationship 
between  the  mean  value  and  standard  deviation  of  the  Richardson  number.  Even 
though  the  above  data  includes  different  synoptic  conditions,  the  relationship 
between  the  turbulence  and  the  Richardson  r.umbe’"  should  be  evaluated  over 
other  types  of  land  surfaces. 


Model  Evaluation 


3  . 

As  an  example  of  the  behavior  of  tne  fl^x  terms  due  to  subgrid  spatial 
correlations,  we  have  iterated  the  three-dimensional,  four  layer,  mesoscale 
model  of  Han  et  a  1 .  (Had)  'see  also  Deardorff  et  al.,  1964)  for  an  idealized 
diurnal  cycle. 

The  entire  mesoscale  model  is  viewed  as  one  grid  box  of  a  larger  scale 
rr.caei  so  that  [4]  is  the  average  value  over  the  entire  mesoscale  model,  ij>  is 
the  local  time-average  evaluated  at  each  grid  point  and  is  the  parameter¬ 
ized  "turbulence" .  The  surface  turbulent  transport  is  again  formulated  with 
the  Louis  relationship  for  the  surface  exchange  coefficient.  The  model 
includes  some  surface  terrain  variations  which  lead  to  nocturnal  drainage  of 
cold  air  for  cloudless  cas-s  of  weak  ambient  flow.  In  the  prototype  numerical 
experiment,  th«  inccminc  solar  radiation  varies  diurnally  leading  to  a  surface 
heat  flux  which  reaches  a  daytime  maximum  of  about  .2  K  m  s_I .  The 
geostrophic  wind  is  specified  to  be  constant  with  a  nominal  speed  of  . 1  m  s_1 
to  simulate  conditions  approaching  free  convection. 

For  the  unstable  daytime  case,  all  c  f  the  spatial  correlation  terms  in 
tne  expression  for  tne  grid-averaged  flux  (4)  are  small  except  for  the 
contribution  due  to  spatial  correlation  between  the  exchange  coefficient  and 
the  surface  wind  speed  'seccnd  term  in  Eg.  4).  This  term  acts  to  reduce  the 
total  grid-averaged  heat  flux  ( [C*7* i<0 ' ,  in  this  case  by  30-40%  (Figure  4a). 

other  words,  wnere  the  wind  speed  is  stronger,  the  instability  tends  to  be 
significantly  less  so  that  t!  exchange  coefficient  is  significantly  smaller. 

However,  the  error  <;ae  to  -  •  -  neglect  of  subgrid  correlations  between 
wind  speed  and  the  exchange  coefficient  is  largely  compensated  by  under¬ 
estimation  of  the  area-averaged  exchange  efficient  (Figure  4b)  which  appears 
in  tne  man  contribution  to  the  ;rui  area-averaged  heat  f 
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Figure  4.  (a)  Various  contributions  to  the  grid-area  surface  heat  flux  (see 

Eq.  4)  for  the  daytime  period.  The  triple  correlation  term  is 
negligible.  (b)  The  total  heat  flux  computed  from  Eq.  4  (open 
circles)  as  compared  to  the  flux  computed  from  area-averaged 
variables  and  area-averaged  exchange  coefficient  (solid  circles) 
and  the  flux  corresponding^ to  the  "model  available"  exchange 
coefficient  computed  from  Ri  (solid  line  with  no  circles; 
essentially  coincides  with  line  with  open  circles  for  this  case). 


£q.  4).  As  a  result,  the  model-estimated  flux  is  close  to  the  true  grid- 
averaged  flux.  Recall  that  the  large  scale  model  can  evaluate  the  exchange 
coefficient  only  in  terms  of  the  Richardson  number  based  on  grid  area-averaged 
variables  (6).  That  is,  the  exchange  coefficient  in  the  large  scale  model 
must  be  computed  as 

Co  =  f (Ri )  (10) 

Ri  =  <g/C0 )( [6 ]-tesfc ])z/[v]2 

where  f  is  the  function  for  the  dependence  of  the  exchange  coefficient  on 
stability.  Compared  to  CH,  the  true  area-average  of  the  exchange 
coefficient  [C^]  is  augmented  by  especiallv  large  values  occurring  at  "hot 
spots"  where  the  instability  is  enhanced  due  to  larger  vertical  gradients  of 
temperature  and/or  weaker  winds.  The  near-cancellation  of  the  important 
spatial  correlation  term  with  errors  due  to  the  underestimation  of  [Cy]  can 
be  shown  to  occur  with  tne  Louis  expression  for  the  conditions  -Ri  >  10  and 
[ 0 1 / 2 ] = [ e jl/2. 

During  the  transitions  between  stable  and  unstable  periods,  the  bulk 
aerodynamic  relationship  in  the  large  scale  model  can  easily  predict  the  wrong 
sign  of  the  grid  area-averaged  heat  flux  as  occurs  in  Figure  5b.  This 
averaging  problem  results  from  the  importance  of  subgrid  correlations  between 
the  exchange  coefficient  and  the  temperature  gradient  (fourth  term  in  Eq.  4). 
Large  upward  heat  flux  in  regions  where  the  vertical  temperature  gradient  is 
unstable  dominates  the  grid  area-averaged  heat  flux  even  though  the  grid 
average  of  the  vertical  temperature  gradient  corresponds  to  stable  stratifica¬ 
tion.  This  "countergradient"  heat  flux  results  from  the  fact  that  the 
exchange  coefficient  is  much  larger  in  the  small  part  of  the  grid  area  which 
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Figure  5.  (a)  Various  contributions  to  the  grid-area  surface  heat  flux  for 

the  nocturnal  period.  (b)  The  total  heat  flux  computed  from  Eq. 
(open  circles)  as  compared  to  the  flux  computed  from  area-averaged 
variables  and  area-averaged  exchange  coefficient  (solid  circles) 
and  the  flux  corresponding. to  the  "model  available"  exchange 
coefficient  computed  from  Ri  (solid  line  with  no  circles).  For 

comparison  with  Fig.  4,  note  the  factor  of  1 03  scale  shift  for  the 
ordinate. 


is  unstably  stratified.  This  particular  averaging  problem  appears  to  be 
rather  short-lived  with  da.urnally  varying  flow.  However,  such  an  effect  could 
exert  a  longer  term  influence  in  situations  which  are  persistently  character¬ 
ized  by  near  neutral  stratification.  This  averaging  problem  is  somewhat 
analogous  to  the  countergradient  heat  flux  resulting  from  time-averaging  the 
vertical  temperature  gradient  in  the  interior  of  the  heated  boundary  layer 
(Deardorff,  1966)  or  time  iteration  which  excludes  diurnal  variations  (Mahrt 
et  al . ,  1986). 

For  the  strongly  stratified,  nocturnal,  situation,  the  relative 
importance  of  the  various  subgrid  correlation  terms  (Figure  5a)  is  quite 
different  from  the  unstable  case.  Now  the  wind  speed  and  the  exchange 
coefficient  are  positively  correlated,  leading  to  enhancement  of  the  downward 
heat  flux.  That  is,  tne  staoility  is  locally  reduced  in  regions  of  strongest 
airflow  which  increases  the  exchange  coefficient. 

This  increase  of  tne  grid-averaged  heat  flux  is  opposed  by  the  negative 
correlation  between  spatial  variations  of  the  exchange  coefficient  and  the 
vertical  gradient  of  potential  temperature.  The  exchange  coefficient  is  small 
where  the  vertical  gradient  of  potential  temperature  is  large.  The  grid- 
averaged  heat  flux  is  also  reduced  by  the  triple  correlation  term  (last  term 
in  Eq.  4)  which  is  physically  quite  complex. 

The  spatial  correlation  terms  in  (4)  sum  to  near  zero  for  the  stable  case 
so  that  the  net  modification  of  the  area-averaaed  heat  flux  due  to  the  subgrid 
correlation  terms  is  small.  That  is,  the  total  flux  is  close  to  that 
predicted  by  the  main  contribution  to  the  area-averaged  heat  flux  (first  term 
in  Eq.  4).  However  this  main  term  is  underestimated  by  models  due  to  the 
fact  that  the  exchange  coefficient  Cg,  which  is  based  on  the  Richardson 
numoer  computed  from  area-averaged  variables,  is  considerably  smaller  than  the 
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true  area-averaged  exchange  coefficient  [CH]  (Figure  5b).  This  under¬ 
estimation  of  the  downward  heat  flux  would  he  expected  to  lead  to  modelled 
surface  temperatures  which  are  too  cold. 

These  results  are  for  the  case  of  strong  radiational  cooling  at  the 
surface.  If  the  geostrophic  wind  speed  is  increased,  subgrid  variations  are 
reduced.  For  winds  on  the  order  of  10  m  s-^ ,  the  averaging  problems  have 
become  negligibly  small  at  least  based  on  the  numerical  experiments  with  the 
mesoscale  model  used  in  this  study.  The  generality  of  these  results  is  not 
known  and  Figures  4-5  must  be  considered  only  as  examples  of  potential 
averaging  problems. 

The  transport  by  subgrid  scale  motions  associated  with  local  topography 
is  not  explicitly  reported  here  because  it  involves  correlations  between 
temperature  and  downslope  flow.  The  bulk  aerodynamic  relationship  describes 
heat  flux  perpendicular  to  the  local  ground.  In  an  absolute  coordinate  system 
where  the  vertical  coordinate  is  parallel  to  the  gravity  vector  and 
independent  of  local  slope,  downslope  currents  lead  to  an  upward  heat  flux. 
Such  a  heat  transport  redistributes  heat  in  response  to  cooling  over  sloped 
terrain.  Since  such  differential  heating  is  not  included  in  numerical  models 
on  subgrid  scales,  the  inclusion  of  such  redistribution  of  heat  would  appear 
to  be  not  appropriate.  The  same  could  be  said  of  the  influence  of  heat 
transport  by  daytime  upslope  currents.  However,  secondary  streets  could  be 
important  particularly  if  the  upslope  currents  initiate  moist  convection. 

The  estimation  of  the  area-averaged  momentum  flux  is  more  difficult 
because  of  the  influence  of  the  terrain-induced  pressure  drag.  This 
pressure  drag  is  often  incorporated  by  enhancing  the  drag  coefficient  or 
surface  roughness  length.  However,  the  practice  of  using  the  suosequently 
enhanced  surface  friction  velocity  in  the  formulations  for  mixing  in  the 
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boundary  layer  is  not  justified  since  the  pressure  drag  of  the  topography  does 
not  translate  into  boundary  layer  turbulence  or  at  least  not  in  a  way  which  is 
described  by  existing  boundary-layer  theory. 
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6. 


Reformulation 


Before  reformulating  t’..e  dependence  of  the  exchange  coefficient  on  the 
Richardson  number,  several  additional  complications  must  be  noted.  In  theory 
the  turbulence  vanishes  as  the  gradient  Richardson  number  exceeds  a  critical 
value.  In  numerical  models,  this  asymptotic  possibility  should  not  be 
invoked,  not  only  because  of  horizontal  averagina  problems  but  also  because 
the  Richardson  number  is  computed  over  a  finite  depth  determined  by  the  model 
resolution.  No  matter  how  large  the  Richardson  number  for  the  grid  layer, 
turbulence  in  actual  atmospheric  flows  can  always  be  generated  over  thinner 
layers  where  the  Richardson  number  is  locally  small.  The  turbulence  over  such 
thin  layers  often  occurs  intermittently  at  cnanging  levels  so  that  flux  over  a 
deeper  layer  is  established  on  a  time  scale  which  is  longer  than  that  of  the 
intermittency . 

In  addition,  momentum  can  be  transported  vertically  by  nonlinear  gravity 
waves  while  significant  vertical  transport  of  heat  may  be  generated  by 
radiational  flux  divergence,  especially  with  strong  surface  inversions.  Both 
of  these  transport  mechanisms  are  generally  neglected  in  large  scale 
boundary-layer  models.  Because  of  these  influences,  and  the  significant 
spatial  averaging  effects  in  the  stable  case,  any  formulation  of  the  surface 
flux  may  suffer  large  errors.  For  this  reason  complicated  schemes  attempting 
to  include  the  details  of  the  influence  of  stability  are  not  justified. 

Both  the  idealized  analysis  _n  Section  3  and  the  specific  calculations  in 
Section  5  indicate  that  for  the  stable  case,  the  exchange  coefficient  relating 
area-averaged  fluxes  to  area-averaged  gradients  should  be  significantly  larger 
than  predicted  by  the  usual  expressions  for  the  exchange  coefficient. 
Furthermore,  the  exchange  coefficient  should  not  vanish  for  large  Richardson 
number  as  is  implied  by  the  observations  reported  in  Section  4  and  as 
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previously  recommended  by  Kondo  et  al .  (1976).  The  case  of  unstable  strati¬ 

fication  does  not  require  systematic  modification  at  least  based  on  the  above 
results,  although  the  flux  due  to  mesoscale  subgrid  motions  could  be 
important . 

A  slower  decrease  of  the  value  of  the  exchange  coefficient  with 
increasing  positive  stability  can  be  most  simply  formulated  as 

C  (Ri)  =  C  exp(-mRi)  (11) 

ri  n(J 

where  CHq  is  the  value  at  neutral  staoility  as  predicted  by  the  Louis 
formulation.  The  idealized  averaging  results  in  Figure  1  suggest  that  m  is  a 
little  greater  than  one  for  unity  standard  deviation  of  Ricnardson  number. 
Considering  tne  dependence  of  the  standard  deviation  on  the  Richardson  number 
(Section  4)  and  additional  influences  which  enhance  the  area-averaged  flux  in 
stable  conditions  (Section  5  and  discussions  above),  m  =  1  appears  to  be  a 
suitable  value. 

Thus,  expression  (11)  for  the  stable  case  and  the  usual  Louis  formulation 
for  the  unstable  case  form  a  tentative  model  of  the  surface  exchange 
coefficient  which  attempts  to  include  the  most  important  qualitative  aspects 
of  grid-area  averaging.  With  present  luck  of  understanding,  (11)  could  be 
applied  to  the  transfer  of  other  quantities  in  addition  to  heat.  The 
performance  of  (11)  in  a  given  large  scale  model  would  presumably  depend  on 
the  details  of  the  model,  especially  the  height  of  the  lowest  model  level.  In 
the  model  of  Troen  and  Mahrt  (1986),  relationship  (11)  applied  to  heat, 
'.omentum,  and  moisture  leads  to  the  expected  enhancement  of  downward  fluxes 
for  the  very  stable  cases  although  realistic  representation  requires  several 
levels  within  the  thin  nocturnal  boundary  layer.  In  models  which  have  been 
indirectly  adjusted  to  compensate  for  the  underestimation  of  downward  heat 
fluxes  and  anomalous  cooling,  the  use  of  (11)  may  not  be  beneficial. 
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7. 


Conclusions 


The  formulation  of  the  subgrid  scale  flux  in  numerical  models  commits 
three  types  of  errors  related  to  the  implied  spatial  averaging  of  the  grid 
area.  First,  the  flux  due  to  subgrid  motions  larger  than  turbulence  scales 
[see  Eq.  (2)]  is  not  included.  We  did  not  examine  this  problem  here  since  it 
is  strongly  dependent  on  situation.  Secondly,  extra  flux  terms  result  from 
the  spatial  averaging  of  the  local  flux-gradient  relationship  [see  Eq.  (4)]. 
These  terms  are  due  to  spatial  correlation  between  the  locally  averaged  varia¬ 
bles  appearing  in  the  flux-gradient  relationship.  Thirdly,  errors  result  from 
necessity  of  relating  the  exchange  coefficient  to  resolved  grid  area-averaged 
variables  instead  of  spatially  averaging  the  local  exchange  coefficient  as 
required  by  the  Reynolds  averaging.  Because  these  errors  can  be  large,  the 
use  of  sophisticated  local  relationships  between  fluxes  and  gradients  dees  not 
appear  to  be  justified  for  use  in  large-scale  numerical  models. 

The  particular,  modeling  results  of  Section  5  indicate  that  the  important 
spatial  correlation  term  for  the  unstable  case  approximately  cancels  errors 
due  to  the  use  of  existing  local  formular.ions  for  the  exchange  coefficient. 

For  the  stable  case,  the  spatial  correlation  terms  approximately  cancel  each 
other  while  the  required  grid-averaged  exchange  coefficient  is  seriously 
underestimated.  This  underestimation  is  due  to  relatively  large  values  of  the 
exchange  coefficient  within  the  parts  of  the  grid  area  where  the  stability  is 
weakest.  This  problem  is  important  because  of  the  strong  nonlinearity  of  the 
relation  between  the  exchange  coefficient  and  stability.  Even  though  the 
absolute  magnitude  of  the  surface  fluxes  is  small  in  the  stable  case,  and 
probably  exerts  little  influence  on  the  overlying  free 'atmosphere,  such  small 
fluxes  become  important  in  the  surface  energy  balance  and  significantly 
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influence  the  surface  air  temperature.  The  convergence  of  downward  turbulent 
flux  of  heat  occurs  over  a  thin  boundary  layer  in  the  strongly  stratified  case 
and  therefore  can  be  locally  significant  even  if  the  flux  magnitude  is  small. 
Then  use  of  the  usual  local  flux-gradient  relationship  and  the  associated 
underestimation  of  the  grid-averaged  downward  heat  flux  will  lead  to 
unrealisticailiy  rapid  surface  cooling. 

A  revised  formulation  (11)  for  the  dependence  of  the  excnange  coefficient 
on  the  Richardson  number  is  constructed  for  the  stable  case.  The  revised 
formulation  is  thought  to  improve  significantly  the  prediction  of  the  grid- 
area  averaged  flux.  However,  any  formulation  for  the  stable  case  remains 
tentative  due  to  the  incomplete  understanding  of  turbulence  with  stable 
conditions  and  due  tc  the  lack  oi  observations  of  spatial  variations  of 
surface  fluxes. 

Observational  verification  for  the  stable  case  is  difficult  since  fluxes 
are  weak  and  computed  flaxes  are  often  seriously  contaminated  by  sampling 
problems.  The  observed  relationship  between  small  scale  vertical  velocity 
variance  and  the  layer  Richardson  number  indicates  three  distinct  regimes  as 
previously  found  in  Kondo  et  al.  (1978),  although  the  values  of  the  Richardson 
number  at  the  transitions  between  regimes  depend  on  the  depth  of  the  layer  for 
the  computation  of  the  Richardson  number.  For  the  weakly  stratified  case 
(Ri  <  1),  the  turoulence  variance  does  not  systematically  vary  with  the 
Richardson  number.  For  the  moderately  stratified  case  (1  <  Ri  <3),  the 
turbulence  strength  decreases  linearly  with  increasing  Richardson  number.  For 
strong  stability  (Ri  >  3),  the  turbulence  is  weak  but  is  not  significantly 
reduced  by  further  increases  of  the  Richardson  number.  Even  with  large  layer 
Richardson  number,  tu.  bulent  tranrport  may  continue  over  thinner  layers,  at 


least  intermittently.  As  a  result,  the  exchange  coefficient  should  not 
totally  vanish  with  large  layer  Richardson  number. 
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CHAPTER  V 


PARAMETERIZATION  OF  THE  STABLE  BOUNDARY  LAYER 


1.  Introduction 

In  this  study,  the  planetary  boundary  layer  (PBL)  model  of  Troen  and 
Mahrt  (1986;  hereafter  referred  to  as  TM86)  is  examined  to  determine  its 
response  to  specific  parameterization  changes.  This  model  is  currently  used 
in  AFGL's  3-D  global  spectral  model  (Brenner  et  al.,  1984;  Yang  et  al., 

1988) .  The  modifications  discussed  below  are  motivated  by  examination  of 
other  model  formulations  as  well  as  analysis  of  aircraft  and  tower  data  in  the 
stable  boundary  layer  ( SBL) .  Any  model  improvements  must  be  weighed  against 
potential  increases  of  computational  costs  and  model  complexity. 

The  parameterization  of  the  very  stable  boundary  layer  is  an  important 
problem.  For  example,  consider  the  forecasting  problems  of  minimum 
temperature  and  pollution  concentrations.  Under  very  stable  conditions  with 
clear  skies  and  calm  winds  at  night,  the  minimum  surface  temperature  can  be 
determined  to  a  large  degree  on  the  solution  of  the  surface  energy  balance 
(Pan  and  Mahrt,  1987)  .  with  inadequate  sensible  heat  transfer,  the  amount  of 
nocturnal  cooling  of  the  surface  can  easily  be  overpredicted.  The  predicted 
formation  of  dew  also  plays  an  important  role  (Oke,  1978)  . 

The  more  statically  stable  the  lower  atmosphere,  the  greater  the 
potential  for  serious  air  pollution  episodes.  The  meteorological  conditions 
which  favor  serious  air  stagnations  are  well-known  to  be  associated  with 
synoptic  scale  high  pressure,  when  large-scale  subsidence  acts  as  a  cap  to 
boundary  layer  growth.  The  presence  of  a  capping  temperature  inversion  on  a 
statically  stable  PBL  restricts  vertical  mixing  of  pollutants.  In  addition, 
horizontal  advection  of  pollutants  away  from  their  source  is  not  likely  near 
the  centers  of  high  pressure  systems,  where  winds  are  weak.  The  mechanism 
which  is  perhaps  most  important  in  this  regard  is  the  strength  of  the 
turbulence  (parameterized  through  the  eddy  diffusivity)  which  can  diffuse  the 
pollutants  upward  from  the  surface.  With  inadequate  vertical  mixing,  the 
model  would  predict  higher  concentration  values  than  might  be  expected. 

The  main  changes  studied  here  involve,  firstly,  the  determination  of  the 
height  of  the  PBL  under  stable  conditions,  where  a  critical  Richardson  number 
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formulation  is  used.  The  second  change  involves  application  of  a  modified 
exchange  coefficient  for  momentum,  heat,  and  moisture,  following  Mahrt  (1987) . 
The  final  change,  and  perhaps  the  most  significant  one,  involves  the 
reformulation  for  the  eddy  diffusivity  in  the  boundary  layer  for  very  stable 
conditions  following  Kendo  et  al . ,  (1978) .  The  motivations  for  considering 

and  ultimately  implementing  these  changes  will  become  apparent  separately  in 
each  of  the  following  sections. 

These  examinations  involve  running  the  one-dimensional  version  of  the 
model  of  TM86,  performing  separate  sensitivity  tests  of  each  of  *-he  new 
formulations.  Finally,  ail  three  changes  are  included  together.  Initial 
>  conditions  for  this  sensitivity  experiment  were  specified  for  a  48  hour  run, 

using  latitude  20°N,  longitude  10°E,  0800  GMT  on  21  June.  The  soil  and 
atmosphere  are  both  taken  to  be  dry,  the  modelled  S'-1’ i  properties  are  those  of 
sand,  the  atmospheric  lapse  rate  is  6°C  km  ^  with  an  initial  surface 
temperature  of  20.7°C;  the  time  step  used  in  the  model  is  180  s  and  the 
vertical  resolution  is  20  m. 

The  suggested  reformulations  are  described  in  detail  in  Section  2. 
Section  3  presents  the  results  of  the  sensitivity  tests  for  the  initial 
conditions  listed  above.  In  Section  4,  other  tests  of  the  1-D  model  are 
discussed,  including  a  simulation  of  conditions  during  the  Wangara  field 
experiment  (Clarke  et  al.,  1971)  and  a  winter  snow  cover  situation  over 
Manitoba  using  a  new  model  package  for  surface  snow  cover.  Finally,  the 
results  are  summarized  in  Section  5. 
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2.  New  Formulations 


a.  PBL  height  and  the  critical  Richardson  number 


The  Richardson  number  is  viewed  as  the  ratio  of  buoyancy  destruction  of 
turbulence  to  its  production  by  shear,  and  in  gradient  form  has  a  theoretical 
critical  value  of  0.25,  according  to  many  studies  over  the  past  fifty  years 
(Panofsky  and  Dutton,  1984)  .  Recently,  however,  the  value  of  the  critical 
Richardson  number  has  again  come  into  question.  Miles  (1987)  (whose  1961 
contribution  had  much  to  do  with  the  instillation  of  the  value  I/4)  gives  an 
historical  accounting  of  the  Richardson  number  and  its  usage  and  finally 
suggests,  based  on  recent  advances  in  nonlinear  hydrodynamic  stability  theory 
(e.g.,  Abarbanel  et  al . ,  1986),  a  value  of  Ric  =  1. 

The  Richardson  number  may  take  different  forms,  but  in  numerical  models 
the  layer  Richardson  number  is  the  only  type  which  can  be  used.  It  is  defined 
as 

_  g  Az  A9 

^  L  2  ’ 

e0[AU] 


with  9q  the  mean  potential  temperature  of  the  layer,  A0  the  change  in 
potential  temperature  across  the  model  layer,  Az  the  thickness  of  the  model 
layer,  and  AU  the  change  in  wind  speed  over  the  layer.  If  the  bottom  layer  is 

the  ground  surface,  AU  becomes  simply  U,  Az  becomes  z,  and  Ri  is  called  the 

D 

bulk  Richardson  number.  Because  of  the  sensitivity  of  Ri  to  the  depth  over 
which  gradients  are  measured,  it  is  generally  believed  that  the  true  critical 
Richardson  number  actually  increases  as  the  depth  increases  (Lyons  et  al., 
1964) .  In  fact,  turbulence  is  often  observed,  at  least  on  an  intermittent 

basis,  for  la^er  Richardson  numbers  much  larger  than  0.25.  A  critical  value 

of  1.0  is  often  used  in  applications  (Brutsaert,  1982),  and  this  value  will  be 
used  in  the  control  experiment. 

In  the  present  model,  the  height  of  the  boundary  layer  for  the  stable 
case  is  determined  from 


Ri  9.,  U2 

c  0 

n  =  - 

g  A9 

°  V 


(2) 
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where  h  is  the  boundary  layer  height,  U  is  the  wind  speed  in  the  model  layer, 
and  A0  is  the  virtual  potential  temperature  change  from  the  top  of  the  layer 
to  a  point  near  the  surface  (see  TMu6»,  and  Ric  is  the  critical  Richardson 
number.  The  determination  of  h  is  an  iterative  procedure,  beginning  at  the 
surface.  At  each  model  level,  the  bulk  Richardson  number  is  calculated  and 
compared  to  the  critical  Richardson  number.  When  Rig  exceeds  Ric,  the 
boundary  layer  height  is  found  by  a  procedure  which  is  illustrated  graphically 
in  rig.  1. 

The  proposed  change  to  TM86  does  not  involve  the  format  for 

determination  of  h  but  rather  the  value  of  the  critical  Richardson  number. 

3ecause  of  the  sensitivitv  of  h  to  Ri  ,  it  is  important  that  a  critical  value 

*  c 

be  chosen  which  simulates  a  realistic  transition  to  suppressed  turbulence  (for 

Ri  >  Ri  ) .  TM86  chose  a  value  of  O.b.  It  is  clear  from  (2)  that  any  change 
c 

in  Ri  will  force  a  direct  response  in  the  predicted  value  of  h.  Observations 
do  support  the  notion  that  turbulence  continues  at  layer  Ri  values  well  above 
1  (e.g.,  Pcrtman  et  ai.,  1962,  Webb,  197C;  Kondo  et  ai.,  1918;  Kunkel  and 
Walters,  1992;  Louis  et  ai.,  1983;  Mahrt,  1987).  In  light  of  these 
observations  and  new  theoretical  developments,  perhaps  a  value  as  large  as  5 
for  a  critical  value  of  the  bulk  Richardson  number  may  not  seem  unreasonable. 
Experiments  will  be  carried  out  for  several  values  of  the  critical  Richardson 
number . 
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b.  S-U.rfa.ce  exchange  coefficients 


As  mentioned  earlier,  Mahrt  (1987)  has  proposed  that  a  modified  exchange 
coefficient  for  heat  be  implemented  in  atmospheric  PBL  models.  Such  a  change 
is  motivated  because  present  models  underestimate  the  strength  of  heat 
transfer  in  the  very  stable  case,  due  to  typical  subgrid-scale  variations  of 
surface  conditions.  Recently,  the  European  Centre  for  Medium  Eange  Heather 
Eorecasts  (ECMWF)  modified  the  parameterization  of  surface  processes  in  their 
operational  model  (Bottger,  1987),  partly  motivated  in  part  by  exaggerated 
nocturnal  cooling  under  clear  skies.  Apparently  PBL  models  often  predict 
surface  temperatures  which  are  too  cold  at  night  due  to  erroneous  elimination 
of  downward  heat  flux  in  the  strongly  stratified  surface  inversion  layer. 

Other  corrections  to  compensate  for  the  lack  of  sensible  heat  transport  under 
statically  stable  conditions  are  sometimes  used,  or  a  stable  layer  may  not  be 
permitted  to  develop,  in  which  case  the  exchange  coefficient  for  the  neutral 
case  is  used,  ensuring  continuation  of  sensible  heat  transport. 

In  TM86,  the  coefficients  of  Louis  (1979)  have  been  used  for  momentum, 
heat,  and  moisture.  Sensitivity  tests  have  been  performed  to  compare  this 
reiationsh.'' p  with  the  one  proposed  by  Marut  for  the  stable  case,  namely 


(3) 


where  is  the  calculated  heat  exchange  coefficient,  is  its  value  for 

neutral  conditions  as  calculated  in  the  Louis  formulation,  and  m  i3  an 
adjustable  parameter  thought  to  be  about  1. 

The  heat  exchange  coefficient  will  be  substantially  larger  for  the  Mahrt 
formulation,  as  can  be  seen  in  comparing  Fig.  2a  to  2b.  Also  as  m  in  (3)  is 
decreased,  the  exchange  coefficient  decreases  more  slowly  with  increasing  bulk 
Richardson  number  (Fig.  2c) .  Results  from  experiments  for  a  few  values  of  the 
adjustable  parameter  m  will  be  presented  in  Section  3.  The  control  run  uses 
the  Louis  formulation  for  . 
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Although  (6)  was  intended  fur  use  in  the  surface  layer,  some  PBL  models  simply 
use  the  same  form  in  the  parameterization  for  turbulence  in  the  outer  layer  as 
well,  with  the  (1  -  z/k)  factor  included  to  allow  for  gradual  decrease  of  K 

m 

as  the  top  of  the  boundary  layer  is  approached. 
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TM 86  show  that  in  their  model,  effectively  reduces  to 


f 

K  =  0.09  L  u*  1 

m  * 

V 


1 . 


(8) 


On  dimensional  grounds,  the  eddy  diffusivity  can  be  seen  a3  the  product  of  a 
velocity  and  length  scale,  so  that 


K 


m 


-  u 


(9) 


Comparing  (9)  to  (8)  would  seem  to  indicate  L  is  one  of  the  relevant  length 
scales  for  turbulent  mixing  above  the  surface  layer  in  TM86.  However,  the 
Monin-Obukhov  length  is  a  valid  length  scale  only  in  the  surface  layer 
according  to  similarity  theory;  the  surface  layer  may  be  only  a  few  meters 
deep  under  very  stable  conditions.  In  some  of  our  1-D  PBL  sensitivity 
experiments,  L  is  about  1  m,  ar i  therefore  the  surface  inversion  layer  extends 
well  above  L.  Brutsaert  (1982)  notes  that  for  Z/T  >  1,  observations  support 

Lt 

use  of  a  constant  value  of  0m,  such  as  that  suggested  by  Kondo  et  al.  (1978), 

<t>  =  6  .  (10) 

m 

Most  recently  Lacser  and  Arya  (1986)  and  Sorbjan  (1987)  nave,  in  observational 

and  modelling  studies,  respectively,  verified  fairly  constant  vertical 

profiles  of  the  nondimensional  shear  functions  0  and  0,  . 

^  m  h 

Since  K  is  inversely  proportional  to  0  ,  a  limitation  of  the  maximum 
n\  m 

value  on  0^  corresponds  to  substantially  larger  dif fusivities,  increased 
mixing  in  the  PBL,  and  a  deeper  surface  inversion  layer,  compared  to  TM86  and 
previous  models .  This  constant  0m  formulation  incorporates  observ-  tions  which 
virtually  suggests  a  third  regime  for  boux.dary  layer  stability  in  the  model, 
namely  the  very  3table  case.  This  tendency  for  three  regimes  of  boundary 
layer  stability  based  on  the  bulk  Richardson  number  has  been  noted  in 
observational  studies  in  the  past,  including  Portman  et  al.,  (1962),  Kondo 
et  al.  (1978),  and  Mahrt  (1987). 

In  the  changed  model  studied  here,  the  third  regime  is  enacted  if 
Z/T  >  1  in  the  outer  layer.  The  new  eddy  diffusivity  is  adopted  here  by 

Li 

modifying  the  nondimensional  shear  function  as  alluded  to  earlier: 
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"-’C 

(1  +  4.^  L-) 

5.7 


L  <  0 


°-r- 


(ii) 


d.  Sur.T.a  ry  or  5°r:5:riv:l  .■  sf.Ji-. 


The  experiments  performed  and  the  parameters  changed  in  the  runs  are 
summarized  in  Table  1;  changes  are  indicated  in  italics.  The  three  changes 
should  show  some  improvement  if  incorporated  into  TM86  separately,  but  we 
might  expect  the  results  to  be  best  if  they  are  ail  coupled  in  one  run.  This 
takes  sense  in  view  of  the  results  of  Suscher  (1987)  and  other  observational 
studies  which  have  shown  that  although  the  large  eddies  in  the  3BL  are 
important  features  of  the  f'.cw,  they  do  not  contribute  much  to  the  overall 
mixing  which,  is  occur:  i..o.  "us  mixing  appears  to  fee  concentrated  at  the 
edges  of  the  eddies.  u  small-scale  feature  cannot  be  adequately  resolved 

in  a  global  model,  so  ccr izcntal  averaging  is  implied. 

The  for...  ..at  ion  of  Mur.rt  (1787)  was  designed  to  represent  the  implied 
horizontal  averaging.  The  increased  surface  exchange  coefficient  leads 
indirect '.y  to  increased  vertical  mixing  through  the  parameter  u>r  as  seen  in 
Fig.  3.  This  figure  is  a  representation  of  only  the  "first-order"  feedbacks 
of  the  parameters  indicated;  others  of  course  are  occurring  which  are  not 
shown . 

An  increase  in  the  surface  temperature  will  result  from  these  last  two 

modifications  described,  and  by  increasing  the  critical  Richardson  number,  we 

can  increase  the  PBL  height.  This  will  produce  mixing  over  a  deeper  boundary 

layer,  perhaps  further  increasing  the  surface  temperature.  Such  an  increase 

of  Ri  is  warranted  also  in  view  of  the  observations,  where  very  small-scale 
c 

mixing  zones  are  seen  to  exist  despite  overall  very  stable  stratification  on 
the  larger  scale.  By  allowing  turbulence  to  continue  in  the  model  boundary 
layer  under  more  stable  conditions  (through  use  of  a  larger  critical 
Richardson  number),  we  are  effectively  adding  another  facet  to  the 
parameterization  of  turbulence  in  the  very  stable  boundary  layer. 
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Table  1 


Sensitivity  Experiments 


Run 

Ric 

ch 

Kh 

CONTROL 

1 . 0 

Louis 

TM8  6 

TM8  6 

0.5 

Louis 

TM8  6 

RIC  =  3 

3.0 

Louis 

TM8  6 

RIC  =  5 

5.0 

Louis 

TM8  6 

MCH1 . 0 

1 . 0 

Mahrt ,  m=*l 

TM8  6 

MCHO .75 

1.0 

Mahrt,  m=0. 75 

TM8  6 

MCHO  .  1 

1.0 

Mahrt,  n-0 . 1 

TM8  6 

KONDOK 

1.0 

Louis 

new  <t^ 

ALL 

5.0 

Mahrt,  m=l 

new 
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Interrelationships  between  several  parameters  in  the  proposed 
SBL  treatment.  Solid  arrows  represent  positive  feedback,  large 
dotted  arrows  indicate  negative  feedback.  The  direction  of  the 
arrows  show  which  variable  is  affected  (at  the  head  of  the 
arrow)  . 


Figure  3  . 


3 .  Results 


a.  Specification  of  PBL  height 

In  the  model  of  TM86,  a  value  of  0.5  foe  the  critical  Richardson  number 
is  used;  the  value  in  the  present  control  run  is  1.0.  Runs  of  the  model  were 
also  made  for  Ric  =  3  and  Ric  =  5.  The  most  pronounced  impact  of  this  change 
is  expected  to  be  found  in  the  PBL  height  from  (2),  this  is  shown  in  Fig.  4. 
The  daytime  boundary  layer  height  depends  very  little  on  the  critical 
Richardson  number;  it  depends  much  more  on  the  temperature  profile.  There  is 
not  a  great  impact  on  other  variables,  such  as  sxm  temperature  (Fig.  5), 
except  that  in  the  TM86  run,  dew  was  predicted,  and  the  additional 
condensational  heating  warmed  the  surface  about  1.5°C. 

The  vertical  temperature  profile  after  46  hr  of  integration  time, 
roughly  the  time  of  minimum  temperature  is  depicted  in  Fig.  6.  This  figure 
will  be  referred  to  in  the  sections  which  follow  involving  other  comparisons. 
The  vertical  temperature  structure  for  the  runs  with  other  critical  Richardson 
numbers  are  qualitatively  similar  to  Rig.  6,  but  there  is  some  tendency  for 
predicted  boundary  layer  height  to  not  match  the  inversion  top,  due  to  the 
sharpness  of  the  initial  potential  temperature  profile.  The  increased 
critical  Richardson  number  does  have  the  desired  effect  of  deepening  the  SBL. 
This  is  important  because  predictions  of  SBL  height  for  Sms1  winds  in  the 
model  of  TM86  seem  unrealistically  shallow. 


PBL  height 


Figure  6. 


TEST  CASE  DRY  INITIAL  STATE  */8  DCG/KU  LAPSE  RATE  CONTROL 
ckd  too rr  ur-  *>  o  long-  iojj  l*vhj-  »  wttul  raft  mo-  «jut-*i  Jtou*-  MUR*  o 
oat*  rus  *i/wt  torch**  orsositi  an  noun  mo*  racorNtNC  • 


Predicted  vertical  temperature  profiles  for  the  control  run 
after  46  hr. 
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This  section  discusses  the  impact  on  the  SBL  formulation  of  changing  the 

formulation  of  the  surface  exchange  coefficient.  Mahrt  (1987)  recommends  a 

value  of  m  in  Eq.  3  or  order  1;  here  we  test  three  values  of  m  (1.0,  0.75,  and 

0.1) .  Because  the  new  falls  off  exponentially  as  -m  Rig,  a  smaller  value 

of  m  will  enhance  the  exchange  coefficient  (fig.  6.2) .  Recall  also  the 

primary  reason  for  this  change;  an  attempt  to  represent  the  implied  horizontal 

grid-area  averaging  of  large  subgrid  scale  variations  of  C^,  which  presumably 

acts  to  increase  the  surface  temperature. 

Inclusion  of  the  modified  exchange  coefficient  increases  the  value  of  C. 

h 

by  a  factor  of  five,  causing  the  skin  temperature  to  increase,  but  only  0 . 1°C 

(fig.  6.7)  The  enhanced  C,  (five  times  the  value  in  the  control  run)  more 

h 

strongly  mixes  the  lowest  layer  so  as  to  reduce  the  20  m  temperature  by  1 . 5°C 
(fig.  6.8) .  Even  with  m=0.1,  the  increase  in  skin  temperature  was  a  modest 
0 . 6°C  compared  to  the  control  run.  Neither  the  PBL  height  nor  the  surface 
energy  balance  terms  showed  any  appreciable  difference  for  the  various  values 
of  m  and  the  control  run. 

Apparently  the  modified  surface  exchange  coefficient  enhances  downward 
mixing  of  heat  when  the  model  boundary  layer  is  stably  stratified,  cooling  the 
atmosphere  but  only  slightly  warming  the  surface  skin  temperature. 
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c.  Enhanced  eddv  diffusivity 

This  section  describes  the  effects  of  the  modified  K^,  which  is  expected 
to  enhance  downward  mixing  of  sensible  heat  at  night,  partially  ameliorating 
the  effects  of  longwave  radiative  heat  loss  from  the  earth's  surface. 

After  46  hr  for  the  control  experiment,  the  boundary  layer  depth  is  only 
40  m  (Fig.  9a),  so  that  there  are  only  three  levels  which  have  non-zero  heat 
flux.  In  addition,  the  heat  flux  at  the  surface  and  lowest  model  level  are 
always  taken  to  be  equal. .  The  effect  of  increased  diffusivity  with  the  Kondo 
formulation  causes  the  surface  heat  flux  to  decrease  from  -2 . 8  W  m  for  the 

_  T 

control  run  to  -10.5  W  m  ^  (Fig.  9b) . 

The  Kondo  formulation  for  eddy  diffusivity  increases  the  predicted 
surface  temperature  by  0 . 6°C  and  increases  the  predicted  air  temperature  by 
several  degrees  (Fig.  10) .  As  a  result,  the  temperature  difference  between 
the  inversion  top  and  the  surface  is  substantially  reduced.  Increased 
magnitudes  of  downward  sensible  heat  flux  lead  to  warming  in  the  lower 
boundary  layer,  and  cooling  at  the  boundary  layer  top  (Fig.  11),  as  expected. 
The  Kondo  modification  apparently  exerts  greater  impact  than  the  other  two 
changes  involving  tunable  parameters. 


110 
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Figure  10.  Time  series  of  sxin  temperature  for  the  enhanced  eddy 

diffusivity  experiment. 


TEST  CASE  DRY  INITIAL  STATE  W/8  DEG /KM  LAPSE  RATE  KONDO  K 
cud  point  ur-  aao  uo*c-  Lprtu-  xs  iwtul  tow  no-  i^or-tuNOUi*  cmw-  o 
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Figure  11. 


Vertical  temperature  profile  after  46  hr  for  the  enhanced  eddy 
diffusivity  experiment. 


d.  Summary  and  synthesis 

The  first  two  revised  formulations  for  the  SBL  parameterization  scheme 
lead  to  the  desired  effects,  an  increased  PBL  height  and  enhanced  surface 
temperature,  although  the  impacts  seen  in  modified  exchange  coefficient 
experiments  were  quite  small.  It  was  decided  that  all  three  formulations 
would  be  incorporated  into  a  final  run;  combined,  these  changes  might  be 
expected  to  work  together  in  a  synergistic  effect. 

By  increasing  the  PBL  height  through  the  modification  of  the  critical 
Richardson  number,  the  Kondo-K  modification  affects  a  much  deeper  boundary 
layer,  significantly  changing  the  vertical  temperature  profile  (Fig.  12).  The 
surface  temperature  is  significantly  changed,  as  well,  increasing  by  nearly 
5°C  over  the  control  run.  The  surface  exchange  coefficient  modification  plays 
a  role  here,  as  well.  The  effect  of  increasing  the  eddy  diffusivity  alone  has 
the  undesirable  effect  of  reducing  the  exchange  coefficient.  Apparently  by 
smoothing  out  the  boundary  layer  temperature  profile  and  momentum  profile 
through  increased  mixing,  u(  will  be  reduced,  since  it  depends  on  the  vertical 

difference  of  wind  speed  between  the  surface  and  lowest  model  layer.  This 
reduction  of  the  surface  exchange  coefficient  does  not  occur  when  all  three 
changes  are  made.  This  effect  is  seen  in  Figure  13,  a  time  series  of  the 
Kondo-K  prediction  of  Ch  and  that  of  the  control  experiment.  The  way  the 
Kondo-K  formulation  works  with  the  other  two  changes  is  reason  enough  for 
recommending  that  all  three  changes  be  implemented  in  further  tests  of  this 
revised  SBL  parameterization  for  the  ID  model.  Their  minimal  impact  on 
daytime  changes  is  important,  as  we  do  ..ot  wish  to  alter  the  character  of  the 
daytime  PBL,  which  is  controlled  by  a  completely  different  set  of  equations 
due  to  the  quite  different  dynamics  involved. 

The  three  modifications  in  concert  lead  to  a  more  realistic  boundary 
layer  depth  for  5  m  s  ^  winds  (172  m  compared  to  20  m;  Fig.  14),  matching  the 
boundary  depth  with  the  inversion  top  during  periods  of  boundary  layer  growth, 
and  warmer  temperatures  due  to  increased  downward  heat  flux. 

The  results  of  all  experiments  are  summarized  in  Table  2  (z^  refers  to 
the  height  of  the  inversion;  T^  is  the  air  temperature  at  the  lowest  model 
level  above  the  surface) .  Of  course,  the  model  initial  conditions  are  highly 
idealized  in  thi3  case.  In  the  next  section,  Wangara  day  33  is  chosen  as  an 
initial  condition  for  one  run,  while  a  winter  snow  cover  simulation  is  carried 
out  in  a  second  simulation.  For  these  comparisons,  only  the  control  run  and 
the  run  with  all  changes  are  shown  here. 
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Figure  12 . 


The  vertical  temperature  profile  after  46  hr  for  the  experiment 
made  using  all  modifications  to  the  SBL  parameterization. 
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Table  2 

Results  of  Changes  to  Stable  Layer  Formulation  for  idealized 

dry  sandy  soil  case 


Run 

h  (m) 

z  i  (m) 

(°C) 

To  1°  C) 

Cont  rol 

40 

38 

2 . 7 

5.0 

TM8  6 

20 

38 

4.3 

9.1 

RIC=3 

142 

38 

2.5 

6.3 

RIC=5 

172 

38 

3.0 

8.4 

Mahrt  C^, 

m=l 

40 

33 

2.8 

3.5 

Mahrt  Ch, 

m=C  .75 

40 

38 

2.9 

3.6 

Mahrt  Ch, 

m=0  . 1 

44 

38 

2.5 

2.9 

Kcndo-K 

42 

38 

3.3 

11.3 

ALL 

191 

181 

8.0 

14.0 

6 


Other  tests 


4  . 


a.  ftoagara  gay  33 

The  initial  vertical  temperature  profile  for  0610  EST  (Eastern  Australia 
Stand  .id  Time)  on  Wangara  day  33  is  shown  in  Fig.  15a.  The  surface 
temperature  is  just  below  freezing  and  a  strong  radiation  inversion  develops 
under  conditions  of  high  pressure,  clear  skies,  and  modest  winds.  The 
modiiied  model  predicts  m_ch  warmer  surface  temperatures  than  the  control  run 
’  Ig.  16) .  The  observed  surface  temperature  at  0604  on  day  35  was  7.3°C.  The 
control  run  predicted  a  surface  temperature  of  2.1°C,  while  the  modified  run 
yielded  7.6°C,  quite  close  to  the  observed  value. 

The  heat  flux  f  r  the  modified  model  run  is  much  larger  than  the  control 
run  (Fig.  17) .  Due  in  part  to  the  increased  heat  flux,  zp  for  the  modified 
run  is  288  m,  while  it  is  only  1.  .  rr  ‘‘or  the  control;  this  is  not  a  result  of 
only  the  enhanced  diffusiv^ty  as  the  Inc..  sed  Ri^  also  plays  a  role.  The 
observed  inversion  top  at  0623  m  ..a,  ur.  day  35  was  400  m  (Fig.  15b);  again 
the  mcd-i-eu  mod:l  results  look  superior.  However,  this  is  tempered  by  the 
observation  tha-  ‘he  strength  of  ths  inversion  (measured  here  in  terms  of  Tz 

i 

-  Tsfc)  for  day  35  was  only  3.4°C;  the  modified  run  had  a  vertical  temperature 
change  of  7.3°C,  while  the  control  run  had  an  inversion  strength  of  13.4°C! 
There  is  certainly  room  for  improvement  and  runs  which  include  prescribed 
vertical  motion  and  perhaps  the  cloud  model  of  Chu  (1986)  are  indicated  in 
future  research.  Nevertheless,  the  comparison  indicates  optimistic  results 
for  the  three  model  changes. 
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height  (m)  height  (m) 


8 


Control 


temperature  (°C) 

(a) 


ALL 


Simulated  Wangara  temperature  profile  for  control  run  after  48 
hr  (a) .  As  in  (a)  except  for  modified  run  (b) . 
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-20 


-3L 


-10 


heat  flux 

(b) 

Simulated  heat  flux  profile  (W  m-^)  after  48  hr  for  the  V 
control  run  (a) .  As  in  (a)  except  for  the  Wangara  modifi 
(b)  . 
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b.  Winter  snow  cover  simulation 

Sensitivity  tests  were  also  performed  for  the  1-D  model  run  with  a 
surface  snow  cover.  The  surface  snow  cover  model  of  Pan  (see  Appendix  C)  is 
included  in  this  run.  The  initial  conditions  for  this  test  are  a  10  cm  snow 
cover,  and  a  1200  GMT  sounding  taken  from  The  Pas,  Manitoba,  during  February 
of  1987.  There  were  northwest  winds  of  about  5  m  s-^-  near  the  surface;  again, 
moisture  effects  are  not  included.  Comparison  to  observations  is  of  course  an 
impossibility  here  because  of  the  neglected  large-scale  forcing. 

A  pronounced  inversion  typical  of  the  arctic  or  well-developed  cP 
(continental  polar)  air  mass  is  present  in  the  initial  temperature  profile 
(Fig.  18) .  The  control  run  is  colder  at  the  surface  than  the  modified  run 

after  46  hr,  as  expected  (Fig.  19) .  The  diagnosed  PBL  height  is  160  m  for  the 

modified  model,  compared  to  41  m  for  the  control  model;  in  each  case  the 

surface  inversion  has  a  deeper  inversion  overlying  it,  a  manifestation  of  the 

original  input  data.  The  vertical  temperature  profile  for  the  modified  case 
(Fig.  20)  is  quite  typical  of  high-latitude  wintertime  inversions  (Dalrymple, 
1966) .  Sensible  heat  transport  seems  to  be  playing  an  important  role  in  the 
maintenance  of  the  inversion  in  this  test. 

Oke  (1978)  presents  a  typical  surface  energy  balance  over  snow  cover  and 
a  main  feature  is  a  negative  sensible  heat  flux.  The  model-simulated  sensible 

-2  -2 

heat  fluxes  are  -3.0  W  m  and  -16.7  W  m  for  the  control  and  modified  runs, 

-2 

respectively;  Oke  reports  values  of  order  -10  W  m  .  Other  than  these 
comparisons,  we  cannot  make  any  definitive  conclusions  either  about  the 
abilities  of  the  snow  cover  model  or  about  the  impact  which  the  modifications 
will  have  on  snow  cover  simulations. 
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5.  Conclusions 


Three  changes  have  been  made  to  the  stable  layer  formulation  of  the 
model  of  TM86  which  significantly  improve  the  performance  of  the  model.  The 
critical  Richardson  number  change  involves  strictly  a  “tunable"  parameter,  in 
that  observations  are  not  definitive.  This  is  probably  the  hardest  to  justify 
based  on  observational  evidence,  but  it  seems  necessary  to  provide  for 
realistically  deep  stable  boundary  layers  with  moderate  wind  speeds.  Strong 
support  exists  for  the  other  two  changes,  however,  enhancements  to  the  surface 
exchange  coefficients  and  boundary  layer  eddy  diffusivity.  In  tests  with  an 
idealized  atmosphere,  the  results  indicate  that  the  three  proposed  changes 
seem  to  work  in  a  synergistic  manner  to  improve  the  prediction  of  SBL  physics, 
including  PBL  depth,  surface  skin  temperature,  and  vertical  temperature 
profile.  Although  improvements  are  noted  for  the  changes  tested  individually, 
only  the  increased  critical  Richardson  number  change  led  by  itself  to  a 
boundary  layer  significantly  different  from  the  control  model. 

The  revised  model  was  also  tested  for  two  very  different  conditions. 

One  was  a  wintertime  snow  cover  situation  from  February  1987  over  the  Canadian 
prairie,  the  other  was  for  the  Wangara  field  experiment.  Although  direct 
comparison  with  observations  is  not  possible  in  either  case  because  this 
simple  model  must  neglect  some  large-scale  processes,  results  of  the  revised 
SBL  model  are  encouraging. 

Once  further  tests  of  these  formulations  are  made,  which  will  be  the 

subject  of  future  research,  other  possible  modifications  to  the  stable  layer 

formulation  might  be  considered.  For  example,  currently  in  TM86  and  in  the 

model  calculations  described  in  the  present  research,  the  ratio  of  the  heat 

exchange  coefficient  (C.  )  to  the  drag  coefficient  (C  )  is  set  to  be  1.35  (this 

n  m 

is  also  true  for  K^/K^) ,  following  Businger  et  al.  (1971).  This  value  (which 

corresponds  to  a  Prandtl  number  Pr  =  1.35-1  =  0.74)  has  been  shown  to  be  valid 
for  near-neutral  conditions,  but  its  value  in  conditions  far  from  neutral 
still  seems  to  be  somewhat  controversial.  For  example,  Kondo  et  al.  (1978) 
suggest  a  ratio  much  smaller  than  one.  One  might  expect  the  logic  behind  this 
is  that  pressure  fluctuations  due  to  gravity  waves,  ubiquitous  in  the  SBL, 
transport  momentum  but  not  heat  (Finnigan  and  Einaudi,  1981,  Caughey,  1982, 
and  Hunt  et  al.,  1985)  .  It  could  be  noted  here  also  that  the  pressure 


perturbations  due  to  the  turbulence  itself  may  lead  directly  to  momentum 
transport  (Schols  and  Wartena,  1986)  .  While  these  observations  are  valid  for 
linear  wave  motions,  the  stably-stratified  PBL  is  thought  to  contain  highly 
nonlinear  motions,  as  the  analysis  of  Ruscher  and  Mahrt  (1987)  and  others  have 
shown;  which  would  transport  heat,  as  well.  Parameterization  of  turbulence  in 
the  free  atmosphere  is  another  area  which  requires  some  research. 

The  results  for  the  snow  cover  simulation  in  particular  should  be  deemed 
very  preliminary  at  this  point.  Comparison  with  observed  data  using  a  more 
realistic  model  (including  vertical  motion,  advection,  and  clouds,  for 
example)  also  needs  to  be  performed.  Still,  the  fact  that  these  simple 
modifications  are  able  to  abate  the  nocturnal  cooling  is  encouraging  and  they 
warrant  further  testing. 
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Appendix  C:  Modelling  the  Snow  Cover 


1.  Introduction 

In  an  effort  to  parameterize  the  boundary  layer  heat  and  moisture 
transport  for  the  AFGL  global  forecast  model  (Brenner  et  al.,  1984),  a 
combined  model  of  the  boundary  layer  and  the  soil  has  been  developed  (Mahrt 
and  Pan,  1984;  Troen  and  Mahrt,  1986;  and  Pan  and  Mahrt,  1987)  .  In  order  to 
implement  the  boundary  layer  package  into  the  global  model,  it  is  necessary  to 
include  the  effects  of  snowcover,  which  can  be  an  important  part  of  the 
boundary  layer-soil  system  (Tuccillo,  1987) . 

Snow  cover  serves  as  the  upper  boundary  of  the  earth's  surface,  thereby 
affecting  the  boundary  layer  as  well  as  the  soil.  Although  snowcover  reduces 
the  available  energy  at  the  surface  because  of  its  high  albedo  for  solar 
radiation  and  high  emissivity  in  the  spectral  range  of  most  terrestrial 
radiation,  its  insulative  properties  exert  the  greatest  influence  on  the  soil 
(Gray  and  Male,  1981)  .  The  thermal  conductivity  of  new  snow  is  roughly  an 
order  of  magnitude  less  than  that  of  most  soils.  As  snow  "ages"  its  thermal 
conductivity  increases,  but  generally  remains  less  than  that  of  most  soil,  and 
it3  albedo  decreases. 

For  the  atmospheric  boundary  layer,  the  presence  of  snow  means  almost 
complete  removal  of  the  soil  heat  flux.  The  nocturnal  cooling  that  i3  usually 
balanced  by  the  soil  heat  flux  (Oke,  1978)  may  lead  to  much  cooler  surface 
temperatures  in  the  presence  of  snow.  Siberia,  northwestern  North  America  and 
Antarctica  are  among  the  regions  where  intense  radiative  cooling  occurs, 
resulting  in  the  formation  of  air  masses  characterized  by  very  low  surface 
temperature  and  strong  surface  inversions  up  to  1  to  2  km  thick. 

In  this  report,  we  will  outline  the  simple  formulation  that  has  been 
developed  to  model  the  snowcover. 

2.  Model 

At  the  present  stage,  the  global  model  does  not  distinguish 
precipitation  as  either  snow  or  rain  as  the  precipitation  prediction  is  based 
on  a  simple  parameterization  method.  We  categorize  fallen  precipitation  as 
snow  when  the  temperature  at  85  kPa  is  below  the  0°C  threshold.  Snow  is  also 
assumed  to  cover  the  grid  box  area  evenly.  The  first  step  in  the  model  is  to 
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make  an  estimate  of  the  snow  heat  flux  G  when  snow  is  present  by  using 
relationship 


s 


where  *cs  is  the  thermal  diffusivity  for  snow,  Ts  is  the  "3kin"  temperature, 

Tsoi’  the  toP_layer  soil  temperature  (in  the  present  model,  the  top  layer 
is  5  cm  thick) ,  and  h^  is  the  depth  of  the  snow  layer  (assumed  to  be  ten  times 
the  water-equivalent  snow  depth) . 

The  thermal  diffusivity  for  snow  depends  on  the  porosity  of  snow  and  can 
vary  from  0.063  W  K  *  m  for  new  snow  with  a  porosity  of  0.95  to  0.71 
W  K  1  m  1  for  packed  snow  with  porosity  of  0.5  (comparable  to  clay) .  Unless 
we  try  to  resolve  the  snow  surface  into  many  layers  and  keep  track  of  the 
"age"  of  each  layer,  we  can  not  possibly  know  the  porosity  of  the  snow  pack. 

In  this  model,  we  choose  the  value  of  0.13  W  K*1  m'1  for  K  which  corresponds 
to  a  porosity  of  0.8.  The  soil  surface  temperature  is  assumed  to  be  the  same 
as  the  top-layer  averaged  soil  temperature.  Thi3  is  supported  by  observations 
that  the  largest  thermal  gradient  below  the  snow  surface  is  near  the  top  of 
the  snow  layer  (Oke,  1978),  due  to  weak  thermal  diffusion  within  the  snow 
layer.  When  snow  falls  over  warm  soil,  the  snow  heat  flux  may  lead  to 
snowmelt . 

The  calculation  of  the  snow  heat  flux  enables  one  to  calculate  the 
potential  evaporation  E^  using  the  surface  energy  balance: 

<l-a)  Si  +  Li  -OT  '4  =  G  +  p  c  C  |V]  (T'-T  )  +  LE  (C2) 

o  P  h  °  P 

where 

Ep  =  PoChIV,(qs(T')-qo)  .  (C3) 

The  terms  on  the  left-hand  side  of  (2)  are  the  downward  short-  and  long-wave 
radiative  flux  and  the  upward  longwave  radiative  flux.  The  terms  on  the 
right-hand  side  are  the  snow  heat  flux,  the  sensible  and  the  latent  heat  flux. 
The  skin  temperature  T '  is  the  temperature  of  the  surface  if  the  snow  surface 

is  evaporating  at  the  potential  rate.  While  the  unit  of  Ep  in  the  surface 

-2  -l 

energy  balance  are  kg  m  s  ,  typical  soil  hydrological  applications  use  the 
units  ms1.  The  conversion  is  accomplished  via  the  density  of  water  (pw  ■» 

10 3  kg  m  3)  . 
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The  snow  evaporates/3ublimates  at  the  following  rate: 


E  - 


h  >  E  At 

s  p 


h  <  E  At 
s  P 


(C4 ) 


When  the  depth  of  the  snow  layer  is  thick,  it  will  evaporate  at  the  potential 
rate  for  an  entire  time  step.  When  the  snow  layer  is  thin  so  that  it  can  not 
maintain  the  potential  rate,  we  assume  the  snow  to  evaporate  evenly  and 
completely  over  the  time  interval  At  . 

Once  the  evaporation  rate  E  is  determined,  the  skin  temperature  is 
calculated  by  solving  the  surface  energy  balance  (2)  again: 

T  -T 

<1-00  Si  +  Li  -CTT  =  D  — - pc  Cl  V|  (T  -T  )  +  LE  .  (C5) 

s  s  ^  r0  p  h  S  O 


If  the  resulting  skin  temperature  is  above  the  melting  point  of  snow 
(Tc  =  273.16  K) ,  the  amount  of  snowmelt  hn  is  calculated  as  follows: 


<1-00  S-L  +  Li  -CTT 

c 


4 


T  -T 

D  — - pc  C.  IV  |  (T  -T  )  +  LE  +  L.h  , 

s  ^  ropr.  cO  lm 

( C  6 ) 


where  L.  is  the  latent  heat  of  fusion. 

i 

When  precipitation  falls  through  the  air  in  the  model,  it  is  assumed  to 
have  the  same  temperature  as  that  of  the  lowest  atmospheric  layer  (this  is 
arbitrarily  assumed  at  the  present  stage) .  Conversion  of  warm  rain  to  ice  or 
snow  may  be  an  important  process  during  warm  front  passages  and  is  included  in 
the  model.  Excess  snow  melt  is  allowed  to  drip  into  the  top  soil  layer  while 
direct  evaporation  from  the  soil  surface  is  inhibited.  Snowcover  therefore 
helps  soil  to  retain  moisture.  Soil  temperature  is  updated  by  assuming  zero 
heat  flux  across  the  snow-soil  interface. 
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VI.  PARAMETERIZATION  OF  SHALLOW  CONVECTION  IN  THE  BOUNDARY  LAYER 


1 .  Introduction 

Shallow  convection  is  thought  to  be  one  of  the  crucial  sub-grid  scale 
mechanisms  for  transporting  moisture  and  heat  from  the  planetary  boundary 
layer  to  the  free  atmosphere.  It  is  also  thought  to  be  an  important  mechanism 
that  provides  a  balance  against  large  scale  sinking  motion  and  maintains  a 
steady  state  well-mixed  subcloud  layer  under  undisturbed  trade  wind 
situations.  By  its  vertical  transfer  of  water  vapor,  it  supports  the  growth 
of  the  deep,  precipitating  cumulus  towers.  Without  this  mechanism  moisture 
will  accumulate  in  the  lower  boundary  layer  and  cut  down  the  latent  heat  and 
sensible  heat  fluxes  from  the  surface. 

Yanai  et  al .  (1973)  introduced  a  budget  approach  to  implicitly  calculate 
the  first  order  sub-grid  scale  turbulent  flux  from  large  scale  observation 
data.  They  defined  Q1  (apparent  heat  source)  to  be  the  heating  due  to  radia¬ 
tion,  the  release  of  latent  heat  by  net  condensation,  and  the  divergence  of 
the  vertical  eddy  transport  of  sensible  heat,  and  Q2  (apparent  moisture  sink) 
to  be  tne  net  condensation  and  the  divergence  of  the  vertical  eddy  transport 
of  moisture.  They  found  the  Q2  profile  in  the  tropical  ocean  to  possess  peaks 
in  the  upper  troposphere  as  well  as  in  the  lower  troposphere.  The  moisture 
sink  in  the  upper  troposphere  is  interpreted  as  due  to  convective 
precipitation.  The  moisture  sink  in  the  ^ower  troposphere  and  relative  min¬ 
imum  in  Q2  between  the  two  peaks  are  interpreted  as  due  to  transport  of  mois¬ 
ture  upward  from  the  lower  troposphere  and  detrainment  of  liquid  water  and 
water  vapor  in  mid-troposphere  by  the  non-precipitating  shallow  clouds.  They 
concluded  that  the  relative  minimum  between  the  two  peaks  which  counteracts 
the  drying  by  environmental  sinking  motion  is  due  to  the  water  vapor  and  li¬ 
quid  water  detrainment  from  the  clouds,  especially  the  shallow  clouds  in  the 
lower  troposphere. 

Nitta  and  Esbensen  (1974)  used  a  similar  approach  to  estimate  the  large- 
scale  heat  and  moisture  budgets  over  the  tropical  Atlantic  Ocean  during 
Phase  3  (22-30  June  1969)  of  the  Barbados  Oceanographic  and  Meteorological 
Experiment  (BOMEX).  For  undisturbed  period  (22-26  June  1969)  with  weak  cumu¬ 
lus  convection,  the  sub-grid  scale  Q1  and  Q2  profiles  both  show  a  minimum  near 
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the  bottom  of  the  trade  wind  inversion  layer.  This  phenomenon  is  interpreted 
as  the  effect  of  moistening  and  cooling  due  to  moisture  detrained  and  re¬ 
evaporated  at  the  top  of  the  shallow  cumulus. 

There  are  some  other  works  that  obtained  similar  results  by  using  Air 
Mass  Transformation  Experiment  (AMTEX)  data  (e.g.  Esbensar.,  1975;  Nitta,  1976; 
Murty,  1976;  Nitta  and  So,  1980).  As  an  example  we  show  Q1  and  Q2  profiles 
derived  from  AMTEX  75  data  in  Fig.  1  which  is  taken  from  Nitta  and  So  (1980). 


Figure  1  Q1 ,  Q2  profile  from  AMTEX75  period  average.  Figure  taken  from  Nitta 
and  So  (1980). 

The  importance  of  the  shallow  convection  is  not  only  demonstrated  by  the 
budget  studies  but  also  in  numerical  modeling  studies.  Tiedtke  (1983)  showed 
that  it  is  necessary  to  parameterize  such  processes  in  the  general  circulation 
model.  He  found  that  with  both  the  Arakawa-Schubert  (1974)  and  the  Kuo  (1965) 
deep  convection  schemes,  moisture  accumulated  in  the  lower  boundary  layer  and 
there  lacks  another  mechanism  to  transport  the  excess  moisture  up  to  the  upper 
air  environment.  By  employing  a  simple  constant  cloud  diffusivity  profile  as 
a  shallow  convection  scheme  in  the  European  Centre  for  Medium-Range  Weather 
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Forecasts  (ECMWF)  operational  model,  Tiedtke  was  able  to  not  only  get  rid  of 
the  excess  moisture  near  the  surface  and  significantly  improve  the  model's 
forecast  but  also  reproduce  the  BOMEX  and  ATEX  (Atlantic  Trade  Wind 
Experiment,  1969)  sounding  data  in  one-dimensional  model  simulations. 

Direct  calculations  of  the  turbulent  heat  and  moisture  flux  (e.g.  LeMone 
and  Pennell,  1976;  Nicholls,  1980)  using  data  from  the  GATE  (GARP  (Global 
Atmospheric  Research  Program)  Atlantic  Tropical  Experiment,  1980),  however, 
have  found  little  evidence  of  the  effect  of  the  shallow  cumulus.  The  study  of 
California  stratus  by  Albrecht  (1985)  using  aircraft  data,  on  the  other  hand, 
does  seem  to  show  an  enhancement  of  the  moisture  flux.  It  is  possible  that 
the  cloud  amount  may  play  a  role  in  the  magnitude  of  the  fluxes  as  the  GATE 
studies  are  for  days  with  little  shallow  convection. 

One  of  the  goals  of  this  work  is  to  derive  a  shallow  convection  scheme 
using  BOMEX,  AMTEX,  GATE,  and  BLX83  (Boundary  Layer  Experiment  1983  in 
Oklahoma,  Stull  and  Eloranta,  1984)  data.  To  be  consistent  with  the  PBL  model 
treatment,  we  seek  to  parameterize  the  effect  of  shallow  clouds  using  a  dif- 
fusivity  formulation  in  which  the  cloud  diffusivity  profile  depends  on  the 
cloud  amount.  The  cloud  amount  forecast  equation  is  derived  from  the  BLX83 
data  and  is  only  relative  humidity  dependent. 

Troen  and  Mahrt  (1986)  developed  a  PBL  column  model  to  parameterize  the 
turbulent  mixing  within  the  boundary  layer  for  a  global  spectral  model  of  the 
Air  Force  Geophysics  Laboratory  (AFGL)  (Brenner  et  al . ,  1984).  Extensive 
tests  of  the  PBL  model  have  been  performed  (Troen  and  Mahrt,  1986;  Pan  and 
Mahrt,  1987)  and  the  results  indicate  that  the  model  is  capable  of  producing 
realistic  simulations  of  the  boundary  layer  under  various  atmospheric  and 
lower  surface  conditions.  The  shallow  convection  scheme  is  developed  for  this 
model.  Results  of  the  sensitive  tests  are  compared  against  the  constant  cloud 
diffusivity  scheme  of  Tiedtke  (1983). 
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2. 


Data 


Data  from  BOMEX,  AMTEX72  and  GATE  were  chosen  for  this  study.  While  the 
data  covered  different  synoptic  situations,  we  chose  to  use  data  during  undis¬ 
turbed  periods  with  a  trade  wind  inversion  capping  the  unstable  boundary  layer 
when  only  shallow  convective  clouds  are  present.  Following  are  some  of  the 
characteristics  of  each  data  set  (see  Table  1  for  summary). 

2.1.  BOMEX 

The  Barbados  Oceanographic  and  Meteorological  Experiment  (BOMEX)  was 
conducted  during  May  and  June  1969  over  the  area  east  of  the  Barbados.  The 
season  was  chosen  to  provide  a  wide  ranqe  of  convective  activities  in  the  ab¬ 
sence  of  well -developed  storms.  The  primary  objective  of  this  experiment  was 
to  determine  the  rate  of  transfer  of  water  vapor,  heat  and  momentum  from  the 
tropical  ocean  to  the  atmosrhere.  This  experiment  is  divided  into  three 
phases.  Only  Phase  3  from  June  22  to  26  is  used  in  this  study  because  of  a 
particular  shortage  of  data  early  on  the  22nd  and  late  on  the  26th.  This  per¬ 
iod  was  marked  by  relatively  undisturbed  trade-wind  weather  with  light  south¬ 
erly  wind  near  the  sea  surface  turning  into  stronger  northerly  wind  aloft. 
Strong  easterly  winds  appeared  in  tne  lower  boundary  layer  and  decreased  slow¬ 
ly  below  the  inversion.  Cloud  cover  was  estimated  from  ATS-2  image- 
enhancement  satellite  pictures  and  is  about  40%  to  50%  during  the  period.  Sea 
surface  temperature  was  28 *C,  the  period-averaged  sensible  heat  flux  was 
15  W  m-2  and  the  latent  heat  flux  was  167  W  m-2  . 

2.2.  AMTEX74 

The  Air-Mass  Transformation  Experiment  in  1974  (AMTEX74)  was  conducted 
during  14  to  28  February  1974  over  the  Kuroshio  region  around  the  Nansei 
Island  of  Japan.  During  this  season  the  cold  air  mass  is  strongly  modified 
when  it  passes  over  the  oceanic  regions  to  the  east  of  the  Asian  continent. 

Due  to  the  large  amount  of  heat  and  moisture  supplied  from  the  sea  surface, 
cumulus  activity  is  strongly  enhanced  as  the  air  travels  farther  over  the 
ocean.  The  primary  objective  of  AMTEX  was  to  clarify  the  role  of  the  cumulus 
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clouds  and  the  boundary  layer  eddies  in  the  air  mass  transformation  process 
and  their  relationship  to  the  development  of  disturbances. 

This  experiment  was  divided  into  three  periods.  Only  the  first  period 
(from  14  to  16  February)  was  used,  during  which  a  local  high  pressure  center 
is  located  to  the  northeast  of  the  AMTEX  area  and  southerly  wind  prevailed  in 
the  lower  layer  below  the  inversion.  During  this  period  the  average  sea  sur¬ 
face  temperature  was  16*C,  cloud  cover  was  45%,  the  surface  latent  heat  flux 
is  252  W  m~  and  sensible  heat  flux  was  42  W  m“* .  Below  the  inversion  south¬ 
erly  wind  speed  increased  with  height  and  turned  into  strong  westerly  wind 
above  the  inversion.  Sinking  motion  was  stronger  than  for  the  other  cases 
used  in  this  study. 

2.3.  GATE 


The  GARP  Atlantic  Tropical  Experiment  (GATE)  was  conducted  during  1974). 
The  experiment  was  divided  into  three  phases  but  only  data  from  Phase  three  is 
used  due  to  shortage  of  data  for  undisturbed  situation?  during  the  first  two 
phases.  During  this  period  the  boundary  layer  vus  capped  by  a  weak  trade  wind 
inversion  and  boundary  layer  itself  was  nearly  neutral  and  only  slightly 
unstable.  Average  cloud  cover  was  only  10%  with  little  cumulus  activity.  The 
Boundary  layer  structure  appeared  well-mixed  for  most  days  so  that  the  modeled 
counter-gradient  effect  becomes  important  for  the  heat  and  moisture  flux 
profiles.  According  to  observation  the  northerly  wind  component  was  stronger 
them  the  easterly  or  westerly  wind  component.  Sea  surface  temperature  was 
about  26°C,  average  latent  heat  flux  was  90  W  m_i  and  sensible  heat  flux  was 
10  W  m~2. 
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Data 

Sensible  Heat 
Flux  (W  m-2 ) 

Latent  Heat 
Flux  (W  m-2 ) 

— 

Cloud 

Cover 

(%) 

SST 

Cc) 

BOMEX 

15 

167 

45 

28 

AMTEX 

42 

252 

45 

16 

GATE 

10 

90 

10 

26 

3. 


The  Model 


3.1.  The  OSU  PBL  Model 

In  the  Oregon  State  University  Planetary  Boundary  Layer  model  (OSU  PBL) 
the  surface  fluxes  are  computed  according  to  the  Louis  (1979)  formula  and  the 
layer  below  50  m  is  considered  to  be  the  constant-flux  layer.  Above  this 
layer,  temperature  (0),  moisture  (q),  and  momentum  (u,v)  fluxes  due  to  the 
boundary  layer  turbulent  mixing  with  environment  are  parameterized  as  dif¬ 
fusive  processes: 


3X  _  3_  f-SX 
3 1  3  z  3  z 


-  Y) 


(1  ) 


where  X  =  [u,v,9,q],  K  is  the  parameterized  coefficient  of  diffusivity  and  y 
is  the  parameterized  counter-gradient  effect  (see  Troen  and  Mahrt,  1986  for 
detail ) . 

Since  this  is  a  1-D  column  model,  there  is  no  horizontal  advection 
process  affecting  the  model  structure  during  the  sensitivity  tests.  We  pre¬ 
scribed  the  initial  vertical  profile  of  horizontal  (u,v)  and  vertical  (w) 
motion,  temperature  (T)  and  mixing  ratio  (q).  Additional  assumptions  are  made 
as  follows: 

1 .  The  large  scale  vertical  motion  field  is  a  constant  profile 
throughout  the  sensitivity  test  period.  The  tumulent  and  convective 
process  cam  only  change  the  wind,  temperature  and  mixing  ratio  fields 
within  the  boundary  layer. 

2 .  Super-saturation  is  removed  by  assuming  that  the  excess  moisture  is 
carried  away  by  the  trade  wind  and  advected  downstream. 

3.  The  sea  surface  temperature  is  assumed  as  a  constant  throughout  the 
sensitivity  test  period. 

4.  Model's  top  is  4  km  which  constrains  the  boundary  layer  top  below 
this  level.  Between  50  m  and  1.05  km,  the  model  resolution  is 
200  m.  Above  1.05  km,  the  model  resolution  is  100  m. 


3.2.  The  OSU  Cloud  Scheme 


A  formula  for  predicting  the  unstable  boundary  layer  shallow  cumulus 
amount  is  derived  from  the  BLX83  aircraft  and  the  AMTEX  sounding  data.  The 


best  parameter  was  found  to  oe  the  relative  humidity  at  the  Lifting 
Condensation  Level  (LCL)  within  the  boundary  layer.  A  least  square  fit  of  a 
quadratic  polynomial  is  used  to  predict  the  cloud  amount  CC: 

CC  =  AO  +  A1  •  RH  +  A2  •  RH2  ,  for  RH  _>  57%  (2) 

and  the  coefficients  calculated  from  the  data  are:  AO  =  -0.7417,  A1  = 

-1.25665  and  A2  =  2.264E-2.  In  Figure  2  we  present  the  observed  relationship 
between  cloud  cover  and  relative  humidity  at  the  cloud  base.  The  solid  line 
in  Figure  2  represents  the  function  in  Eq.  (2). 

Slingo  (1980)  also  developed  a  cloud  parameterization  scheme  derived 
using  GATE  data  for  use  in  the  British  Meteorological  Office's  tropical  model. 
Slingo's  formulae  for  low  cloud  cover  (Fig.  2)  with  and  without  inversion  are: 

(a)  For  (d9/dp)  <  -0.07(K/mb),  we  assume  the  existence  of  inversion 

min  — 

and  use  the  formula 

CC  =  -16.67 (d9 /dp )  .  +  6 (RH-80 )2 /400  -  1.167  ,  (3) 

mm 

where 


r  1  RH  >  80% 
0  RH  <  80% 


to  estimate  the  cloud  amount. 

(b)  For  (d0/dp)  >  -0.07(K/mb) ,  we  calculate  the  cloud  amount  as 

min 

r  (RH-80)2/400,  for  RH  <  80% 

*  0,  for  RH  <  80%. 


For  the  case  with  inversion  we  use  d9/dp  =  -0.1(K/mb).  For  comparison,  we 
take  the  BOMEX  undisturbed  period  when  d9/dp  is  -0.11(K/mb)  and  cloud  cover  is 
about  50%.  The  relative  humidity  just  below  the  inversion  is  at  91%  and  at 
the  LCL  it  is  84%.  In  this  case  the  Slingo  formula  predicts  a  cloud  cover  of 
70%  while  our  formula  predicts  55%  cloud  cover.  It  cam  be  seen  that  the 
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Slingo  formula  will  overestimate  the  cloud  amount  in  the  presence  of  an 
inversion  wnen  the  relative  humidity  falls  below  80%. 

Because  the  observed  moisture  fluxes  from  either  budget  studies  (AMTEX 
and  BOMEX)  or  aircraft  measurements  (California  stratocunulus  experiment 
reported  in  Albrecht,  1985)  represent  the  combined  fluxes  from  both  the  boun¬ 
dary  layer  turbulence  and  the  shallow  convection,  we  need  to  consider  both 
mechanisms  to  model  the  observed  fluxes .  Within  the  framework  of  the  PBL 
parameterization  scheme,  it  is  decided  that  a  diffusion  type  parameterization 
scheme  will  be  used  to  model  the  moisture  flux  due  to  shallow  convection.  We 
therefore  seek  to  model  the  moisture  flux  w'q'  as 

-f'S.  *  KcK|f  -  r)  (5) 

where  is  the  coefficient  of  diffusivity  due  to  boundary  layer  turbulence 
and  K__  is  the  coefficient  of  diffusivity  due  to  the  shallow  cumulus  within  the 
boundary  layer.  Furthermore,  the  boundary  layer  turbulent  flux  is  calculated 
using  the  Troen  and  Mahrt  (1986)  method.  The  combined  diffusivity  (K^  +  Kc) 
(Fig.  3)  is  first:  calculated  for  each  data  set  from  profiles  of  w'q'  and  q. 

We  observe  a  peak  in  the  lower  boundary  layer  diffusivity  profile.  Troen  and 
Mahrt  (1986)  demonstrated  that  the  profile  of  K  with  a  maximum  at  z/h  =  1/3 
fits  the  profile  derived  from  large-eddy  simulation  experiments  (Wyngaard  and 
Brost,  1983).  By  assuming  that  the  diffusivity  coefficient  for  shallow  con¬ 
vection  (X^)  to  be  small  near  ground,  we  estimate  the  diffusivity  coefficient 
for  the  boundary  layer  (K^)  using  the  profile  of  calculated  K  (=  +  K^) 

below  .5  z/h.  When  a  maximum  K  value  is  observed  near  the  level  z/h  »  1/3,  it 

is  used  to  obtain  a  vertical  profile  of  the  K,  .  We  can  further  derive  a  dif- 

n 

fusivity  profile  for  the  cloud  (Kc>  after  subtracting  from  K  (see 
Fig.  4).  We  found  that  this  profile  can  be  approximately  fitted  by  a  Gaussian 
Distribution  function  (see  Fig.  5a,  b)  if  we  know  the  peak  value  of  this 
function  and  the  standard  deviation.  Based  on  the  observation  that  the  ratio 
of  the  height  of  the  Lifting  Condensation  Level  (LCL)  and  the  boundary  layer 
height  (h)  is  about  0.3  and  that  K  approximately  vanishes  above  z/h  = 

1.2,  we  thus  parameterize  the  profile  of  cloud  diffusivity  as 

K  =  K  exp[-(ZH  -  Center )Z/52 ],  0.3  <  ZH  <  1.2  (6) 

c  max  —  — 


> 
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California  2 


California  1 


Figure  3.  Vertical  distribution  of  (Kc  +  K^)  calculated  from  the  BOMEX, 
AMTEX  and  California  coast  stratocumulus  experiment  data. 


z/h 


(b) 


Figure  5. 


Examples  of  the  best  fit  K  profile  calculated  from  observational 
data.  (a)  For  the  BOMEX,  AMTEX74,  (b)  for  the  AMTEX75. 
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where 


62  =  {1.0  -  Center)2  , 

Center  =  Z  /h  , 
max 

K  is  the  peak  value,  ZH  is  2/h,  Z  is  the  height  where  the  maximum  value 
nicix  max 

took  place,  and  Center  is  about  0.75  for  the  data  we  used  here. 

There  are  several  parameters  (stability,  Richardson  number,  humidity, 

etc.)  that  have  been  tested  to  predict  the  magnitude  of  the  variable  K  , 

max 

but  only  the  cloud  cover  is  significantly  related  to  K  .  The  relationship 

max  r 

between  them  can  be  fitted  by  a  least-square  linear  equation  (Fig.  6)  as 
below: 

K  =  A1  •  CC  +  A2 ,  for  CC  >  40%  ,  (7) 

max  — 

where  A1  =  1.872  and  A2  =  -63.59.  We  further  assume  a  linear  relation  for 
cloud  cover  less  than  40%  as  follows: 

K'  =  A3  •  CC  .  (8) 

max 

where  A3  =  0.3555.  Sensitivity  experiments  reveal  that  simulation  results  are 

insensitive  to  whether  the  K'  term  is  included  or  not.  It  is  probably 

max 

because  the  term  i3  usually  quite  small. 

In  order  to  apply  the  shallow  convection  parameterization  in  real  data 

situations  it  is  necessary  that  Eq.  (6)  can  represent  a  variety  situations 

when  K  may  exist  at  different  levels  within  the  boundary  layer.  Since 
max 

we  are  only  interested  in  the  physical  process  within  the  boundary  layer,  we 
assume  that  cloud  diffusivity  vanished  when  ZH  is  greater  than  1 .2  or  smaller 
than  0.3.  The  ratio  between  the  LCL  height  and  the  boundary  layer  height  is 
close  to  the  0.3  in  all  cases  we  examined  and  the  cloud  top  is  usually  about 
the  same  height  as  the  boundary  layer.  Thus,  we  can  rewrite  5  as  a  function 
of  the  LCL  height  by  assuming  that  the  distance  between  zmax  and  LCL  versus 
the  distance  between  the  boundary  top  and  LCL  is  a  constant  ratio.  We  then 
arrive  at  the  equation: 

Center  =  z  /h  =  9/14  +  5/14  •  z,  .  ,  (9) 

max  lcl 

where  is  height  of  the  LCL. 
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Thera  are  no  similar  treatments  to  the  heat  and  momentum  fluxes  partly 
because  it  is  not  clear  at  the  present  how  the  cloud  redistributes  heat  and 
momentum.  Observational  evidence  seems  to  indicate  that  the  effect  of  shallow 
cumulus  on  heat  and  momentum  budget  is  weak.  We  are  therefore  neglecting  them 
for  the  time  being.  Preliminary  results  show  that  if  heat  and  momentum  are 
similarly  treated  as  moisture,  boundary  layer  will  grow  and  the  inversion  will 
be  wiped  out  within  a  very  short  time  period.  Further  studies  are  necessary 
to  examine  the  effect  of  the  shallow  cumulus  in  the  heat  and  momentum  budget. 


3.3.  The  EC.MWF  PBr,  and  Shallow  Convection  Parameterization  Scheme 


The  European  Centre  for  Medium  Range  Weather  Forecast  (ECMWF)  planetary 
boundary  layer  parameterization  scheme  (Louis,  1979)  does  not  have  an  explicit 
treatment  for  unstable,  stable  or  neutral  boundary  layer.  It  uses  a 
Richardson  number  to  determine  each  grid  level's  instability : 


„ .  g  •  dz  •  d0 

Ri  =  - - x - - 

0 ( dv ) 2 


(10) 


One  then  applies  the  following  formula  to  compute  the  diffusion  coefficient  K 


and 


*h 


where 


1  - 


*  l2  | dv/dz j  F(Ri) 
b  •  Ri 


(11) 


1  +  c  iRil 1/2 


,  unstable  (Ri  <  0) 


F(Ri) 


- stable  and  neutral  (Ri  >  0), 

il  +  b'  •  Ri)2  ~ 


l  = 


k  •  z 


1  + 


k  •  z 


and 


l2  «  b  •  \  (z  +  dz)/z1//3  -  l] 


c  =  C* 


3/2 


.1/2 


dz 


3/2 


The  parameter  F(Ri)  is  a  stability  function,  k  is  the  Von  Karman  constant  (0.4 
in  this  study),  l  is  the  mixing  length,  A  is  the  asymptotic  mixing  length  and 
i3  adjustable  (currently  being  chosen  as  100  m),  b  =  9.4,  b'  =  4.7  and  C*  is 
7.4  for  momentum  and  5.3  for  heat  and  moisture. 
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The  shallow  convection  scheme  used  in  Teidtke  (1983)  is  quite  simple.  A 
constant  diffusivity  coefficient  (25  m2  s-1  )  is  assumed  throughout  the  entire 
cloud  deck  for  momentum,  heat  and  moisture.  The  cloud  base  is  assumed  to  be 
the  condensation  level  for  surface  air  and  the  cloud  top  is  the  level  of  non¬ 
buoyancy,  but  not  higher  than  750  mb.  In  the  OSU  PBL  model  the  cloud  top  is 
determined  where  the  relative  humidity  falls  below  57%  (Eq.  2). 
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4. 


Results 


4.1.  Real  Data  Simulations  with  the  OSU  Scheme 

In  this  section  we  will  present  the  results  of  a  1 -D  model  simulation  of 
the  undisturbed  trade  wind  situations  and  compare  the  simulated  atmospheric 
profiles  with  and  without  the  proposed  shallow  cumulus  parameterization  scheme 
to  study  the  impact  of  the  new  scheme.  The  maintenance  of  the  trade  wind  in¬ 
version  is  generally  thought  to  be  due  to  a  balance  between  the  large-scale 
sinking  motion  of  the  subtropical  high-pressure  system  and  the  turbulent  mix¬ 
ing  within  the  boundary  layer.  To  simulate  this  balance,  it  was  necessary  to 
include  the  observed  vertical  motion  in  the  simple  1-D  model.  Without  large- 
scale  advection,  however,  the  model  predicted  temperature  field  near  the  sur¬ 
face  will  in  time  approach  the  temperature  at  the  sea-surface.  The  virtual 
heat  flux  from  the  ocean  and  the  boundary  layer  turbulence  will  both  be 
reduced.  These  results  will  eventually  lead  to  the  collapse  of  the  boundary 
layer.  Given  this  inherent  limitation  of  the  1-D  model,  we  will  only  present 
short-range  (1-2  days)  mode  simulations. 

For  a  given  initial  profile  of  wind,  temperature  and  moisture,  we  keep 
the  model  top  level  wind  field  and  the  entire  vertical  motion  field  fixed  with 
time  and  allow  bovindary  layer  mixing  to  modify  the  rest  of  the  parameters. 
Because  the.  observed  atmospheric  profiles  are  usually  stably  stratified,  the 
simulated  boundary  layer  will  grow  vertically  as  turbulent  mixing  starts.  We 
will  refer  to  the  initial  state  as  the  time  0  3tate  and  the  predicted  state  by 
the  number  of  hours  from  the  initial  time.  In  addition  to  the  temperature  and 
dew-point  profiles,  we  will  also  display  the  saturation  point  (SP)  profile 
(Betts  and  Miller,  1984).  When  the  SP  profile  is  close  to  a  straight  line,  as 
is  often  observed,  the  atmosphere  is  well -mixed.  In  the  1-D  model,  turbulent 
and  shallow  cumulus  mixing  are  the  only  mechanisms  that  can  modify  the  SP 
profile.  We  have  found  the  SP  profile  to  be  a  good  tool  to  examine  the  effect 
of  the  parameterization  schemes . 

Figures  7  and  8  show  the  model  predicted  atmospheric  profiles  using  BOMEX 
June  22  to  June  26  average  sounding  as  the  initial  state  after  one  12-hour  run 
with  and  without  the  shallow  convection.  By  comparing  these  figures  we  find 
that,  after  12  hours,  both  cases  maintain  a  well-detined  inversion  at  the  top 
of  the  boundary  layer.  This  is  because  sinking  motion  is  warming  the  top  of 
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Figure  8.  Temperature,  dew  point  and  saturation  point  profile  for  the  BOMEX 
case  without  the  OSU  shallow  convection  scheme  at  hour  12. 
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the  inversion  and  cooling  the  upper  part  of  the  boundary  layer  (discussion 
about  the  cooling  will  be  given  in  section  4.2). 

The  temperature  profiles  in  both  Figs.  7  and  8  are  very  close  to 
dry-adiabatic.  The  dew-point  profiles  fall  closely  on  a  constant  mixing-ratio 
line.  The  boundary  layer  coefficient  of  diffusivity  formulation  slightly 
overestimates  the  mixing  near  the  top  of  the  boundary  layer  and  results  in  a 
mixing  ratio  profile  that  increases  slightly  with  height  above  500  m.  Above 
the  boundary  layer,  the  large  scale  sinking  motion  is  adiabatic  and  does  not 
alter  the  SP  character  so  that  the  environmental  SP  profiles  are  maintained. 
The  effect  of  the  shallow  cumulus  parameterization  is  quite  small  at  this 
time,  when  shallow  cumulus  is  parameterized,  the  enhanced  diffusion  will 
transport  more  moisture  into  the  upper  boundary  layer  and,  in  fact,  will 
transport  some  moisture  above  the  inversion.  This  is  designed  to  simulate  the 
penetrating  tops  of  some  of  the  cumulus  that  can  exist  above  the  inversion. 

Above  1.5  km,  the  model  structures  are  identical  as  they  should  be. 

The  boundary  layer  is  significantly  deeper  for  the  case  without  shallow  cumu¬ 
lus  (1.3  km)  than  with  cumulus  (1.0  km).  This  is  because  the  cloud  moisture 
diffusivity  above  the  boundary  layer  top  (h)  mixes  the  dry  stable  air  into  the 
boundary  layer  and  creates  a  deeper  transition  layer  below  the  inversion  top. 
Between  the  inversion  base  and  1.5  km,  the  model  with  shallow  cumulus  (Fig.  7) 
is  actually  moistened  while  the  model  without  shallow  cumulus  (Fig.  8) 
actually  dries  out  because  of  imperfect  vertical  advection  (an  upstream  scheme 
is  applied  and  the  constant  sinking  motion  has  a  maximum  below  1.5  km). 

The  boundary  layer  in  the  1.0  -  1.3  km  interval  is  saturated  for  the  run 
without  shallow  cumulus  (Fig.  8).  This  result  demonstrates  that  the  param¬ 
eterized  boundary  layer  mixing  alone  is  capable  of  transporting  significant 
amounts  of  moisture  into  the  upper  boundary  layer.  This  also  can  be  see  from 
the  model  diagnostic  Q2  profiles  shown  in  Fig.  9  where  the  moisture  flux 
convergence  beyond  one  hour  near  t’  e  top  of  the  boundary  layer  is  greater  in 
the  case  without  shallow  convection  than  with  shallow  convection.  This  is 
because  in  the  latter  case  the  dryer  air  is  mixed  down  from  above  the  top  of 
the  boundary  layer  and  moisture  is  transported  out  of  the  boundary  layer;  this 
can  prevent  moisture  from  accumulating  in  the  upper  boundary  layer.  For  this 
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Figure  9.  Model  diagnostic  Q2  profiles  for  the  BOMEX  case  at  hours  4,  8,  12. 

(a)  With  the  OSU  shallow  convection  scheme,  and  (b)  without  the  OSU 
shallow  convection  scheme. 


Figure  9.  (Continued) 


run  the  predicted  cloud  cover  for  both  cases  is  reasonable  (the  averaged 
amount  is  about  50%  compared  to  45%  for  the  observation). 

4.2.  Comparison  with  Observations 

Judging  by  the  moisture  flux  profiles  deduced  from  the  observed  Q2 
profiles,  shallow  convection  indeed  has  had  an  influence  on  the  environment . 

In  the  1-D  model,  we  also  notice  differences  in  the  SP  profiles  when  shallow 
cumulus  is  parameterized  indicating  changes  in  the  environment.  The  model  re¬ 
sults  further  demonstrate  that  the  often  observed  minimum  in  the  Q2  profiles 
(Fig.  10)  near  the  base  of  the  inversion  is  due  to  the  convergence  of  boundary 
layer  turbulent  flux  as  well  as  due  to  shallow  convection.  In  cases  when  the 
cloud  amount  is  small,  boundary  layer  turbulent  mixing  can  explain  all  or  most 
of  the  observed  Q2  profile.  When  shallow  cumulus  is  included  in  the  1-D 
model,  the  predicted  minimum  in  Q2  is  nearly  twice  the  observed  minimum  in 
magnitude.  This  result  can  be  explained  by  the  fact  that  the  model-predicted 
moisture  profiles  become  saturated  after  a  few  hours.  As  the  moisture  gradi¬ 
ent  across  the  inversion  increases,  the  parameterized  turbulent  moisture  flux 
also  increases.  In  a  recent  effort  to  include  the  OSU  boundary  layer  model  in 
a  global  spectral  primitive-equation  model  (to  be  reported  elsewhere),  the 
boundary  layer  structure  over  the  ocean  is  closely  monitored  and  is  rarely 
found  to  be  saturated.  It  is  obvious  that  mechanisms  such  as  horizontal 
advection,  convective  heating  and  large-scale  precipitation  are  also  important 
in  determining  the  boundary  layer  moisture  profiles. 

In  Table  2  we  list  the  surface  latent  heat  flux  and  virtual  heat  flux  at 
hour  12  from  the  model  calculation  for  the  three  data  sets  (unit  is  W  m“2 ) . 
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Figure  1 Oa 


Vertical  profiles  of  large-scale  apparent  heat  source  Q1 , 
apparent  moisture  sink  Q2,  radiation  heating  QR  for  AMTEX  during 
14-16  Feb.  1974.  Thin  dashed  line  denotes  the  inversion  base. 


Figure  10b  Vertical  profiles  of  apparent  heat  source  Q1 ,  apparent  moisture 
sink  Q2,  and  radiation  heating  for  the  BOMEX  undisturbed  period. 
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Table  2.  Summary  of  the  surface  latent  heat  and  virtual  heat  flux  simulated 
with  the  OSU  scheme.  (Note:  *  represents  without  shallow 
convection . ) 


Virtual  heat  flux 
(W  m-2 ) 

Latent  heat  flux 
(W  m"2) 

BOMEX 

13  (9)* 

194  (188)* 

AMTEX 

155  (153)* 

675  (648)* 

GATE 

3* 

55* 

Comparing  the  model  results  with  the  observations  (Table  1 ) ,  we  find  that 
in  general  the  model  generates  higher  latent  heat  flux  and  lower  sensible  heat 
flux.  This  is  because  saturated  sounding  results  in  an  overestimation  of  the 
cloud  cover  and  enhancement  of  the  cloud  d^ffusivity  as  well  as  the  surface 
moisture  flux.  The  virtual  heat  flux  is  slightly  underestimated.  Because  the 
potential  temperature  increases  with  height,  enhancing  the  heat  transport  due 
to  shallow  cumulus  would  lead  to  a  warmer  temperature  in  the  lower  boundary 
layer.  This  would  further  reduce  the  sensible  heat  transport  lending  further 
support  for  not  parameterizing  the  heat  transport.  The  increase  in  the 
virtual  heat  flux  with  shallow  cumulus  (Table  2)  is  due  to  an  increase  in  the 
turbulent  moisture  transport  which  results  in  a  slightly  drier  lower  boundary 
layer  and  a  higher  evaporation  rate  from  the  ocean. 

In  the  AMTEX  simulation  both  fluxes  are  overestimated  at  hour  12  and  in 
the  GATE  case  both  fluxes  are  underestimated.  However,  as  discussed  in  the 
previous  section,  for  the  AMTEX  data  set,  the  model  takes  12  hours  to  adjust 
itself  toward  a  steady  state.  After  that  the  surface  moisture  flux  became 
302  W  m~  and  sensible  heat  flux  became  37  W  m  .  For  the  GATE  data  set,  the 
model  does  not  include  the  cumulus  effect  because  the  LCL  is  above  the  boun¬ 
dary  layer. 

4.3.  Comparison  with  the  ECMWF  Scheme 

Here  we  would  like  to  compare  the  boundary  layer  parameterization  scheme 
and  the  shallow  convection  scheme  with  the  ECMWF  schemes.  Unlike  our  boundary 
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layer  parameterization  scheme,  the  Louis  formula  used  in  the  ECMWF  model  de¬ 
pends  only  on  the  Richardson  number  [Eq.  (10)].  It  is  not  necessary  for  the 
diffusivity  coefficients  or  [Eq.  (11)]  to  go  to  zero  near  h  (top  of 
the  boundary  layer),  although  the  normal  stability  of  the  atmosphere  is  such 
that  the  and  become  very  small  above  the  boundary  layer. 

A  recent  ECMWF  report  (Louis  et  al . ,  1981)  shows  two  modifications  of  the 
Louis  formula  which  may  increase  the  boundary  layer  turbulent  mixing  because, 
in  day-to-day  diagnostic  of  the  global  operational  model,  it  became  apparent 
that  the  model  suffers  from  a  lack  of  boundary  layer  mixing.  The  first  modi¬ 
fication  includes  a  larger  asymptotic  mixing  length  (X)  of  300  m  to  increase 
the  and  K^.  The  second  modification  not  only  slightly  changes  the 
formulation  but  also  increases  the  asymptotic  mixing  length  to  400  m  for  heat 
and  moisture.  The  modification  is  in  response  to  artificial  cooling  of  the 
stratosphere  generated  by  the  model  and  strong  weakening  of  the  jet  stream  in 
the  model  integrations.  Essentially,  the  increased  asymptotic  mixing  length 
for  the  heat  and  moisture  will  enhance  the  turbulent  mixing  processes  and  re¬ 
distribute  the  heat  and  moisture  into  higher  region.  Furthermore,  Tiedtke 
(1983)  developed  a  simple  shallow  convection  scheme  for  ECMWF  in  response  to 
deficiencies  found  in  global  integrations  which  is  attributed  to  a  lack  of 
shallow  convection. 

In  order  to  illustrate  the  difference  between  the  OSU  and  the  ECMWF 
boundary  layer  parameterization  scheme,  the  coefficients  of  diffusivity  using 
both  models  are  presented  in  Fig.  11  (first  modification  of  the  Louis  scheme 
is  used  for  ECMWF).  The  initial  state  is  the  BOMEX'case  and  the  ECMWF  scheme 
is  used  to  predict  for  24  hours.  At  hour  24  the  coefficient  of  diffusivity 
based  on  our  scheme  is  also  calculated.  It  can  be  seen  that  the  OSU  paramet¬ 
erization  scheme  produces  a  profile  of  that  is  about  an  order  of  magnitude 
greater  than  that  from  the  ECMWF  scheme.  In  addition,  the  from  the  ECMWF 
scheme  vanishes  above  .7h  (the  be  dary  layer  height  determination  is  based  on 
the  OSU  scheme).  This  is  a  result  of  the  chosen  asymptotic  mixing  length  that 

oroduces  a  maximum  for  K,  around  300  ra.  Stronger  mixing  in  the  OSU  scheme 

h 

leads  to  warmer  and  drier  surface  air.  This  induces  larger  surface  latent 
heat  and  sensible  heat  fluxes. 

Here  we  implement  both  Louis'  and  Tiedtke 's  schemes  in  our  boundary  layer 
routine.  Figures  12  and  13  are  the  result  from  BOMEX  initial  state.  Each 
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Figure  11.  Vertical  profile  of  the  moisture  diffusion  coefficient  for  the  OSU 
(dash  line)  and  the  ECMWF  (solid  line)  schemes. 


of  these  runs  uses  asymptotic  mixing  length  of  300  m.  Due  to  the  weak 
turbulent  mixing  process  as  was  shown  above,  it  takes  more  than  24  hours  for 
both  ECMWF  schemes  to  adjust  the  boundary  layer  structure  from  the  initial 
state  to  well-mixed.  Compared  to  our  result  (less  than  one  hour  can  result  in 
a  well-mixed  lower  boundary  layer)  the  turbulent  mixing  for  ECMWF  scheme  is 
relatively  weak  and  slow.  Hence  we  display  model  structures  at  hour  36. 

Table  3  lists  the  results  at  hour  36  from  the  ECMWF  scheme. 

Table  3.  Summary  of  the  surface  latent  heat  and  virtual  heat  flux  simulated 
with  the  ECMWF  scheme.  (Note:  *  represents  without  shallow 
convection . ) 


Data 

Virtual  heat  flux 
(W  m-2  ) 

Latent  heat  flux 
(W  m-2 ) 

BOMEX 

1.5  (0.9)* 

32  (18)* 

GATE 

0.45  (0.42)* 

17  (13)* 

The  characteristics  of  the  simulated  BOMEX  model  structure  are:  1)  the 
predicted  boundary  layer  height  is  lower,  2)  a  steady  state  structure  can  not 
be  maintained  for  very  long,  and,  3)  the  upper  boundary  layer  is  not  saturated 
with  or  without  the  shallow  convection.  With  the  shallow  convection  a  well- 
defined  inversion  layer  is  obtained.  Folding  of  the  SP's  profile  above  the 
inversion  indicates  the  transition  from  the  moist  surface  layer  to  the  dryer 
mixed  layer,  but  this  is  not  very  clear  for  the  cases  without  the  shallow  con¬ 
vection  due  to  the  weaker  turbulent  mixing. 

Closer  examination  of  the  potential  temperature  (6)  and  mixing  ratio  (q) 
profile  can  reveal  more  detail  of  the  boundary  layer  model  structure.  For 
example,  for  the  BOMEX  case  at  hour  36  a  deep  mixed  layer  with  a  top  at  800  m 
is  created  with  shallow  convection;  thi3  is  roughly  twice  the  depth  simulated 
without  the  shallow  convection.  Unlike  our  deeper  and  totally  well-mixed 
boundary  layer  model  structure  at  hour  12  (see  Fig.  8),  the  ECMWF  result 
within  the  boundary  layer  shows  a  thick  moist  surface  layer  below  250  m  which 
explains  their  smaller  simulated  surface  virtual  heat  flux  and  latent  heat 
flux  (Table  3). 


Down-wind  of  the  trade  wind  cumuli  is  the  typical  area  where  we  cam  see 

the  deep  cumulus.  It  is  very  important  to  the  development  of  the  down-wind 

deep  convection  to  transport  moisture  from  the  surface  to  the  upper  boundary 
layer.  With  the  ECMWF  scheme,  however,  moisture  is  trapped  in  the  lower  boun¬ 
dary  layer  and  cannot  be  transported  to  the  higher  region  without  the  shallow 

convection.  A  comparison  of  Pig.  13  to  Fig.  8  shows  that  the  dew  point  de¬ 

creases  rapidly  above  500  m  for  the  ECMWF  scheme,  while  our  scheme  simulates  a 
1.1  km  depth  well-mixed  layer  with  a  nearly  constant  dew  point.  This  may  be 
the  reason  why  the  ECMWF  global  model  needs  the  shallow  convection  to  set  up  a 
deeper  mixed-layer  and  enhance  the  surface  turbulent  fluxes.  In  Fig.  8  our 
result  indicates  that  the  shallow  convection  only  slightly  increases  the 
moisture  above  1.3  km  and  the  profile  below  that  level  is  not  changed.  But 
from  the  potential  temperature  point  of  view,  the  shallow  convection  creates  a 
deeper  transition  layer  above  the  boundary  layer  and  the  top  of  the  boundary 
layer  becomes  lower. 

The  primary  difference  between  the  schemes  can  be  seen  in  the  model 
diagnostic  Q2  profiles  at  hour  36.  For  the  BOMEX  case  (Fig.  14),  the  Q2  pro¬ 
file  without  shallow  convection  shows  a  large  moisture  flux  convergence  near 
the  surface  and  a  weaker  one  at  the  top  of  the  mixed  layer  (400  m).  When 
shallow  convection  is  included  the  Q2  profile  shows  only  one  minimum  at  the 
mixed  layer  top  (800  m).  This  means  the  modeled  shallow  cumulus  is  working  to 
prevent  moisture  from  accumulating  near  the  surface  and  to  transport  it  from 
the  lower  boundary  layer  to  the  mixed  layer  top.  This  also  enhances  the  sur¬ 
face  sensible  and  latent  heat  transport.  Comparing  the  case  with  shallow  con¬ 
vection  to  our  result  (which  essentially  is  without  shallow  cumulus),  we  find 
that  both  have  a  minimum  point  at  250  m  and  the  magnitude  is  slightly  larger 
for  the  OSU  model. 

For  the  BOMEX  case  with  shallow  convection,  the  ECMWF  scheme  has  a 
reasonable  Q2  minimum  compared  to  the  observation.  The  boundary  layer  top  (h) 
is  much  lower  and  the  minimum  point  of  the  Q2  profile  is  found  at  the  top  of 
the  mixed  layer  instead  of  the  boundary  layer  top.  However,  the  simulated 
thick  moist  surface  layer  is  not  observed  in  the  real  data.  In  contrast,  our 
simulation  has  a  deeper  and  well-mixed  boundary  layer  (which  is  close  to  the 
observation)  and  a  slightly  overestimated  moisture  flux  at  the  top  of  the 
boundary  layer. 


Q2(c/day) 


Figure  14.  Predicted  model  diagnostic  Q2  profiles  for  the  BOMEX  case  at  hours 
4,  8,  12.  (a)  With  the  ECMWF  shallow  convection  scheme,  and 

(b)  without  the  ECMWF  shallow  convection  scheme. 


From  the  above  discussion  one  can  see  that  the  role  of  the  shallow 
convection  is  crucial  to  the  ECMWF  model  because  the  boundary  layer  mixing  is 
too  weak.  By  including  the  shallow  convection  in  the  model,  the  turbulent 
mixing  is  still  not  capable  of  reproducing  the  state  that  is  close  to  the  real 
data.  Tftese  are  due  to  the  constraint  of  the  asymptotic  mixing  length  that 
prevents  moisture  and  heat  from  being  mixed  higher  into  the  atmosphere. 


5. 


Conclusion  and  Discussion 


The  ECMWF  global  prediction  model  is  considered  by  many  the  state-of- 
the-art  weather  prediction  model  at  the  present.  It  is  the  result  of  recent 
model  diagnosis  from  the  ECMWF  that  led  to  our  interest  in  modeling  the  shal¬ 
low  cumulus.  However,  results  based  on  our  boundary  layer  parameterization 
scheme  are  quite  different  from  the  ECMWF  experience.  We  find  that  our  boun¬ 
dary  layer  parameterization  can  mix  moisture  into  the  upper  troposphere  and 
create  the  observed  Q2  profile  while  the  ECMWF  scheme  cannot.  The  results 
indicate  that  the  primary  mechanism  that  transports  moisture  away  from  the 
lower  atmosphere  is  the  boundary  layer  turbulent  flux.  The  boundary  layer 
turbulent  mixing  alone  is  capable  of  maintaining  an  apparent  moisture  source 
near  the  inversion.  While  the  sensible  heat  flux  over  the  ocean  becomes  quite 
small  after  a  few  hours,  the  virtual  heat  flux  remains  positive  and  the  boun¬ 
dary  layer  remains  in  the  unstable  regime. 

Due  to  the  constraint  of  the  asymptotic  mixing  length,  we  have  found  the 
ECMWF  boundary  layer  mixing  is  restricted  to  the  lowest  500  m  of  the  boundary 
layer.  Even  with  recent  modifications,  the  boundary  mixing  is  still  extremely 
small  above  500  m.  It  is  possible  that  the  use  of  the  diffusivity  coefficient 
formula  derived  for  a  neutral  boundary  {Blackadar,  1962)  in  unstable  situa¬ 
tions  may  be  inappropriate. 

The  effect  of  the  shallow  convection  scheme  in  our  1-D  model  is  to 
enhance  the  boundary  layer  turbulent  mixing  and  the  surface  turbulent  fluxes 
and  to  reduce  moisture  flux  convergence  near  the  top  of  the  boundary  layer  by 
mixing  the  air  within  the  boundary  layer  with  the  free  atmosphere.  For  the 
ECMWF  model,  the  shallow  convection  scheme  significantly  improves  the  model 
results  by  enhancing  the  surface  turbulent  fluxes  as  well  as  the  moisture  flux 
convergence. 

It  is  necessary  to  derive  a  proper  boundary  layer  parameterization  scheme 
before  one  tries  to  study  the  importance  of  the  trade  wind  shallow  cumulus. 
Without  doing  so  one  can  be  very  easily  misled  by  the  results  and  would  have  a 
wrong  picture  of  the  problem. 

This  study  is  primarily  focused  on  the  trade  wind  shallow  cumulus.  In 
the  future  we  would  like  to  see  the  effect  of  the  boundary  layer  and  shallow 
convection  scheme  over  other  areas  (e.g.,  over  land  situation)  and  would  also 
like  to  focus  on  its  impact  on  the  3-D  model. 
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CHAPTER  VII:  CONCLUSIONS 


The  work  under  the  present  contract  has  concentrated  on 
improvements  of  generic  problems  in  boundary-layer  modelling  such  as 
transport  induced  by  shallow  cumulus,  subgrid  variations  of  surface 
fluxes,  the  unique  behavior  of  transport  within  the  very  stable 
boundary  layer  and  interaction  between  soil  hydrology, 
evapotranspiration  and  boundary-layer  development.  This  work  is 
described  in  some  detail  in  the  current  report  and,  except  for  the 

[  shallow  cumulus  effort  (see  discusion  below) ,  is  either  already 

> 

published  or  soon  will  be  submitted  for  publication  in  major  journals. 

The  present  work  has  attempted  to  develop  physically-motivated 
formulations  which  avoid  practical,  but  ad  hoc,  corrections.  We 
therefore  anticipate,  without  proof,  that  our  new  formulations  for 
shallow  cumulus  transport,  the  surface  exchange  coefficient,  and 
transport  within  the  very  stable  boundary-layer  are  not  as  model 
sensitive  as  ad  .ioc  corrections  and  would  therefore  be  more  robust 
with  respect  to  any  future  major  changes  in  the  rest  of  the  model. 

Considering  our  emphasis  on  somewhat  independent  improvements  of 
various  aspects  of  the  boundary-layer  ;  ackage,  the  overall  package 
seems  to  be  surprisingly  compatible.  However  conclusive  statements 
cannot  be  made.  In  fact,  incompatibility  problems  are  almost  sure  to 
arise  as  the  total  boundary-layer  package  is  further  studied  within  the 
Air  Force  Global  Spectral  Model  under  a  variety  of  different 
meteorlogical  conditions.  Some  possible  problem  areas  might  include 
the  interaction  between  the  snow  physics  and  the  boundary-layer  model 
during  long  periods  of  strong  surface  radiative  cooling,  or,  the 
interaction  between  the  representation  of  boundary-layer  cloud 
transport,  surface  evaporation  and  boundary-layer  growth  with  a  variety 
of  cloud  situations.  The  "real"  improvement  of  the  present 
developments  within  the  global  model  must  be  further  re-evaluated  when 
more  realistic  respresentations  of  the  global  distribution  of  surface 
properties  are  adopted. 


More  specifically,  the  impact  of  the  boundary  layer  on  the  global 
circulation  must  be  studied  in  detail  for  different  geographic  regions 
and  different  synoptic  situations.  This  extensive  homework  is 
necessary  because  the  interaction  between  the  boundary  layer  and  the 
free  atmosphere  involves  complex  nonlinear  coupling.  As  a  result  the 
impact  of  the  boundary  layer  mouel  is  not  always  what  it  seems;  the 
intuitive  zero  order  effects  are  often  exceeded  by  secondary  effects 
not  previously  anticipated.  The  analysis  of  boundary  layer-free  flow 
interactions  must  also  address  the  appropriateness  of  the  soil  and 
vegetation  specification  and  the  initialization  of  soil  water. 
Considerable  future  effort  will  be  devoted  to  four-dimensional  data 
assimilation  which,  includes  updated  soii  moisture  and  temperature. 

A  major  shortcoming  of  the  present  model  is  the  absence  of 
t irbulenc  or  subgrid  scale  transport  above  the  boundary  layer.  In  the 
real  atmcspher-.-,  the  turbulence  is  sometimes  stronger  above  the 
boundary  layer.  This  is  particularly  frequent  in  the  case  of  the  very 
stable  nocturnal  boundary  layer  where  turbulence  is  sometimes 
maintained  in  weakly  stratified  layer  corresponding  to  the  mixed 

layer  .; :  the  previous  daytime  period.  Other  examples  include  the 
f-'rmaticn  of  clear  air  turbulence  and  strong  turbulence  in  certain 
flews  over  complex  terrain. 

Therefore  a  formulation  is  requi-ed  which  allows  for  local  shear- 
generation  of  turbulence  and  topographi-  .lly  induced  transport  of 
momentum  and  ether  quantities  associated  with  nonlinear  gravity  waves. 
An  important  aspect  of  this  problem  is  the  effective  Prandtl  number  of 
the  formulated  transport.  Nonlinear  gravity  waves  and  turbulence  in 
stratified  flow  usually  lead  to  larger  values  of  the  momentum 
diffusivity  as  compared  to  the  diffusivity  for  heat  and  other 
quantities.  This  difference  is  due  to  the  generation  of  momentum 
transport  by  pressure  effects.  Obviously  this  problem  is  also 
important  in  the  very  stable  boundary  la^er. 

The  boundary  layer  cloud  formulitior.  requires  father  examination. 
This  formulation  plays  a  crucial  role  in  the  surface  energy  balance, 


the  overall  boundary-layer  development  and  the  influence  of  the 
boundary  layer  on  the  general  circulation.  The  problem  is  difficult 
because  acceptably  simple  approaches  necessarily  emphasize  either 
cumulus  or  stratus  type  boundary  layer  clouds,  and  the  behavior  of 
radiative  and  evaporative  cooling  at  the  cloud  top  occur  on  scales  that 
cannot  be  resolved  in  the  present  model.  Future  work  will  concentrate 
on  modification  of  the  enhancement  of  the  diffusivity  due  to  shallow 
cloud  convection.  The  present  formulation  overestimates  transport  as 
the  total  cloud  cover  (including  inactive  clouds)  approaches  100 
percent . 

The  most  challenging  stage  of  the  future  research  will  be 
evaluation  of  boundary  layer-general  circulation  interactions  with  the 
improved  shallow  cloud  formulation  and  the  inclusion  of  free 
tropospheric  transport  by  turbulence,  gravity  waves  and  other  subgrid 
scale  processes.  This  work,  must  be  carried  out  jointly  with  the 
Atmospheric  Prediction  Branch  of  the  Air  Force  Geophysics  Laboratory. 

We  have  recently  been  granted  a  new  contract  from  AFGL  to  continue  our 
research  in  these  directions. 
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ABSTRACT 


4 


A  two-iayer  model  of  soil  hydrology  and  thermodynamics  is  combined 
with  a  one-dimensional  model  of  the  planetary  boundary  layer  to  study 
various  interactions  between  evolution  of  the  boundary  layer  and  soil 
moisture  transport .  Boundary  layer  moistening  through  surface 
evaporation  reduces  the  potential  and  actual  surface  evaporation  as  well 
as  the  boundary-layer  growth.  With  more  advanced  stages  of  soil  drying, 
the  restricted  surface  evaporation  allows  greater  sensible  heat  flux 
which  enhances  boundary- layer  growth  and  entrainment  drying. 

Special  individual  cases  are  studied  where  the  wind  speed  is 
strong,  solar  radra-icn  is  reduced,  transpiration  is  important,  the  soil 
is  thin,  or  the  soil  is  covered  with  organic  debris. 
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ABSTRACT 


This  study  examines  the  inadequacies  of  formulations  for  surface 
fluxes  for  use  in  numerical  models  of  atmospheric  flow.  The  difficulty 
is  that  numerical  models  imply  spatial  averaging  over  each  grid  area. 
Existing  formulations  are  based  on  the  relationship  between  local  fluxes 
and  local  gradients  and  appear  to  poorly  describe  the  relationship 
between  the  grid-averaged  flux  and  the  grid-averaged  gradient.  For 
example,  area-averaging  the  bulk  aerodynamic  relationship  reveals 
additional  spatial  correlation  terms  and  a  complex  relationship  between 
the  grid-averaged  exchange  coefficient  and  the  stability  based  on  "model 
available"  grid-averaged  variables. 

This  problem  is  studied  by  assuming  idealized  spatial 
distributions  of  the  Richardson  number  over  a  grid  area.  Some 
perspective  is  provided  by  consulting  observed  spatial  distributions  of 
the  layer  Richardson  number  at  the  surface.  Various  contributions  to 
the  area-averaged  surface  flux  are  studied  by  employing  a  small-scale 
numerical  model  as  a  grid  box  of  a  larger-scale  numerical  model.  Based 
on  these  analyses,  a  new  formulation  is  proposed  for  relating  the  area- 
averaged  flux  to  the  area-averaged  gradient.  However,  this  expression 
cannot  be  seriously  tested  with  existing  observations. 
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ABSTRACT 


The  very  stable  boundary  layer  is  a  reaicr.  of  the  atmosphere 
typified  by  large  vertical  gradients  of  temperature  and  momentum. 
Analysis  of  very  stable  atmospheric  flows  is  complicated  by  the  presence 
of  nonlinear  interactions  among  gravity  waves,  shear-driven  overturning 
circulation,  two-dimensional  vortical  modes  and  intermittent  turbulence 
ir.  various  stages  of  development.  This  study  examines  the  horizontal 
structure  of  a  very  stable  atmospheric  boundary  layer,  using  data 
obtained  primarily  from  terrain-following  aircraft  flights  over  central 
Oklahoma . 

Several  diagnostic  procedures  are  applied  to  the  aircraft  data, 
including  classical  and  rotary  spectral  analysis,  principal  component 
analysis,  and  structure  functions.  Coherent  structures  with  sharp 
boundaries  are  examined  with  a  new  conditional  sampling  technique  which 
requires  little  a  priori  specification  of  sampling  criteria.  Because 
the  flows  involve  sharp  boundaries,  spectral  techniques  do  not  provide 
as  much  useful  information  as  other  more  localized  procedures.  The 
edges  cf  the  coherent  structures  are  regions  of  significant  vertical 
heat  transport,  a  feature  not  often  emphasized  in  studies  of  gravity 
waves  and  vortical  modes  in  the  stable  boundary  layer. 

The  presence  of  significant  turbulence  even  for  large  stability 
has  implications  for  modelling  of  the  very  stable  boundary  layer. 
Forecasts  of  minimum  temperature,  boundary  layer  height,  inversion 
characteristics,  and  pollutant  dispersal  are  all  significantly  affected 
by  turbulent  mixing.  Many  models  of  the  stable  boundary  layer 
artificially  arrest  the  mixing  under  stable  conditions,  resulting  in, 
for  example,  overestimates  of  nocturnal  cooling.  A  new  parameterization 
of  the  stable  boundary  layer  is  studied  here  by  incorporating  it  into  an 
existing  model  of  the  planetary  boundary  layer.  The  model  is  then  run 
with  one-dimensional  sensitivity  tests  for  an  idealized  atmosphere  and 
with  data  from  Wangara  day  33.  A  simulation  over  snow  cover  is  also 
examined.  The  tests  substantiate  the  role  of  vertical  mixing  in 
ameliorating  nocturnal  cooling,  validating  the  parameterization  changes. 
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ABSTRACT 


A  snailow  convection  scheme  is  derived  from  several  data  sets 
(BOMEX,  GATE,  AMTEX,  BLX83)  AND  DEVELOPED  FOR  THE  OSU  1-D  boundary  layer 
model.  Results  of  the  model  structure  and  characteristics  of  the 
saturation  point  (SP)  profile  are  compared  against  the  constant  cloud 
diffusivity  scheme  of  Tiedtke  (1983)  and  the  ECMWF  boundary  layer 
parameterization  scheme. 

The  results  indicate  that  the  primary  mechanism  that  transports 
moisture  away  from  the  lower  boundary  layer  is  the  boundary  layer 
turbulent  flux  ar.d  that  the  boundary  turbulent  mixing  alone  is  capable 
of  maintaining  an  apparent  moisture  source  near  the  inversion.  While 
the  sensible  heat  flux  over  the  ocean  becomes  quite  small  after  a  few 
hours  of  model  simulation,  the  virtual  heat  flux  remains  positive  and 
the  boundary  layer  remains  in  the  unstable  regime. 
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